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POST  STABILIZATION  IONIZATION  LFA'KL  PKKDICTIONS 
Volume  III  of  the  Calendar  Year  1975  Annual  Report  to  the 
Defense  Nuclear  Atjenfy 

Section  1 
INTRODUCTION 

Widespread  ionization  of  the  lower  atmosphere  by  emissions  from 
radioactive  debris  clouds  formed  in  the  aftermath  of  high  altltvide 
nuclear  explosions  is  potentially  a source  of  interference  to  satellite 
communication  systems.  Debris  clouds,  advected  by  mesospheric  winds, 
have  been  shown  by  Zalesak  and  Coffey  1075^  to  spread  over  distances 
of  hundreds  of  kilometers  from  the  burst  point  within  a few  hours  after 
detonation  in  addition  to  undergoing  dramatic  changes  in  shape.  Typi- 
cally, a cloud  develops  from  an  initial  spherical  shape  into  a highly 
elongated  configuration  under  the  influence  of  short-scale  length 
vertical  shears.  The  predictive  capability  for  communications  inter- 
ference is  very  much  a function  of  the  accuracy  of  the  mesospheric 
wind  field  data  available. 

The  wind  fields  used  by  Zalesak  and  Coffey  (1975^  3re  given  by 
Groves  (1969),CIRA  (1972),  based  on  observational  data  acquired  over 
many  years.  The  mean  deviations  associated  with  these  data  Indicated 
the  presence  of  inconsistencies  that  were  confirmed  by  a subsequent, 
more  detailed  analysis.  This  analysis  of  the  observational  data,  which 
appears  in  section  2 of  this  volume,  provides  the  impetus  for  the  NRL 
program  to  develop  theoretical  models  of  the  mesospheric  circulation 
which  could  be  used  to  substantially  improve  the  wind  data  reliability. 

Initial  studies  of  the  mesospheric  circulation  were  undertaken 

as  a part  of  last  year's  program  with  the  assembly  of  a linear  model 

for  the  mean  circulation  and  both  analytic  and  numerical  simulation 
Note:  Manuscript  submitted  April  7,  1977. 
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of  solar  tidal  phenomena  which  represent  substantial  diurnal  and  semi- 
diurnal perturbations  superposed  on  the  mean  winds.  Further  progress 
in  each  of  these  areas  is  reported  in  this  volume. 

The  mean  circulation  of  the  mesosphere  is  driven  primarily  by 
the  differential  heating  caused  by  the  absorption  of  solar  radiation 
by  ozone.  The  parametization  of  radiative  heating  is  comparable  in 
importance  to  the  parametization  of  Rayleigh  friction  and  radiative 
cooling  which  were  optimized  in  the  initial  NRL  linear  model (Baker  and 
Strobelj  1975).  The  radiative  heating  function  taken  from  the  work 
of  Leovy  (1964)  has  been  replaced  by  the  more  accurate  one  described 
in  section  2 of  this  volume.  Results  from  the  updated  model  containing 
the  new  radiative  heating  parametization  are  also  presented. 

The  response  of  the  atmosphere  to  diurnal  and  semidiurnal  heating 
has  been  studied  at  NRL  with  a three-dimensional  numerical  model 
initially  reported  last  year.  The  model,  which  uses  pressure  as  a 
vertical  coordinate,  has  been  improved  by  insertion  of  a prognostic 
equation  to  determine  the  lower  boundary  condition  at  every  time  step 
in  contrast  to  the  earlier  versions 's  assumption  that  the  substantive 
derivative  of  the  pressure  was  zero  at  the  lower  boundary.  Grid  re- 
solution has  also  been  increased  to  about  0.5  km  from  the  4 km 
resolution  of  the  earlier  version.  This  permits  the  more  accurate 
simulation  of  propagating  waves  under  the  general  prescription  that 
eight  to  ten  grid  points  per  wavelength  are  necessary. 

A fast  running  code  suitable  for  systems  applications  has  been 
prepared  which  follows  the  movement  of  a debris  cloud  using  input 
wind  data  and  calculates  the  ion  density  distribution  attributable 
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to  beta  emission  from  the  radioactive  debris.  The  betas  are  guided 
by  the  geomagnetic  field  lines  into  regions  conjugate  to  the  cloud. 

A description  of  this  code  along  with  a discussion  of  the  accuracy  of 
the  beta  deposition  treatment  appears  in  section  h. 

The  work  performed  at  NRL  during  the  past  year  is  summarized  in 
the  remainder  of  this  volume.  Much  of  this  work  has  been  previously 
reported  in  NRL  Memo  Reports^  at  several  symposia,  in  technical  journals, 
and  at  DNA  sponsored  meetings.  The  principal  contributors  in  each 
technical  area  are  listed  as  authors  of  the  section  describing  that 
work. 
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Section  2 


A CRITICAL  ANALYSIS  OF  CLIMATOLOGICAL  WIND  DATA 

USED  IN  ITiE  FORECAST  OF  RADIOACTIVE  DEBRIS  CLOUD  MOVEMENT 

M.  R.  SCHOEBERL 

It  is  of  considerable  importance  to  communication  system 
performance  ;ELFj  VLF,  HF ) in  a nuclear  environment  to  be  able  to  predict 
debris  cloud  transport  in  the  mesosphere.  Zalesak  and  Coffey  (1975) 
have  shown  that  transport  can  profoundly  alter  the  location  of 
radioactive  debris  clouds.  The  subsequent  beta  decay  within  the  cloud 
results  in  long  lasting^  widespread  ionization  that  can  severely  affect 
the  performance  of  communication  systems.  Accurate  prediction  of 
system  performance  depends  critically  on  our  ability  to  describe  the 
movement  of  debris  patches  by  mesospheric  wind  systems. 

The  purpose  of  this  report  is  to  assess  the  present  method  by 
which  the  spread  of  radioactive  debris  clouds  in  the  upper  atmosphere 
is  determined  and  to  suggest  guidelines  for  future  research.  The 
current  technique  used  to  forecast  debris  cloud  advection  uses  a 
Lagrangian  computer  code  with  model  wind  fields  ( Zalesak  and  Coffey^ 

1975)»  The  wind  fields  are  given  by  Groves  (I969),  CIRA  (1972). 

These  wind  fields^  being  based  upon  many  years  of  observational  data, 
represent  the  best  statistics  available  at  the  time  the  study  was 
undertaken.  The  results  produced  by  Zalesak  and  Coffey  clearly  hinge 
upon  the  accuracy  of  the  wind  field  data.  It  is  thus  important  to 
examine  the  data  carefully  and  ask  if  these  wind  fields  actually  give 
an  accurate  representation  of  the  upper  atmospheric  motions  relevant  to 
determining  debris  cloud  advection.  Groves'  data  are  an  "average"  or 
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climatological  representation  of  the  wind  field  from  which  some  spatial 
and  temporal  fluctuations  have  been  removed  by  the  averaging  process. 

If  these  fluctuations  are  quite  small,  then  the  averaged  data  may  be  used 
to  accurately  forecast  the  transport  of  trace  constituents.  On  the  other 
hand,  if  the  fluctuations  are  large,  the  actual  wind  field  may  rarely 
resemble  the  climatological  wind  field  and  resulting  debris  cloud 
forecasts  based  upon  the  latter  will  be  reliable  only  in  a climatological 
sense. 

In  this  report  the  wind  fields  given  by  Groves  shall  be 
examined  with  the  following  criteria.  First,  we  can,  to  some  extent, 
quantitatively  assess  variability  within  the  wind  field  data  by  examining 
the  s d deviation  of  the  climatological  average  published  by  Groves 

(r 

becond,  we  examine  the  consistency  of  the  data  with  theoretical 
models  of  upper  atmospheric  dynamics.  While  inconsistency  between 
empirical  models  and  theoretical  models  does  not  necessarily  imply 
unreliable  data,  consistency  allows  us  to  use  theoretical  models 
where  empirical  data  may  be  lacking  or  difficult  to  obtain.  From  this 
viewpoint  we  can  determine  if  the  wind  structure  of  the  upper  atmos- 
phere as  computed  from  theoretically  postulated  heat  and  momentum 
sources  bears  any  resemblance  to  the  empirical  wind  structure  given 
by  Groves.  Or  conversely,  we  can  compute  the  implied  heat  and 
momentum  sources  required  to  maintain  Groves'  wind  model  and  compare 
with  known  sources.  Both  aspects  of  this  problem  will  be  discussed. 

Third,  the  zonal  wind  model  of  Groves  is  tested  for  stability 


to  small  wave  perturbations.  Instability  probably  implies  the  presence 
of  large  scale  eddy  mixing  which  could  greatly  affect  the  transport 
of  radioactive  debris. 

Within  the  body  of  this  report  the  data  and  its  observed 
variability  of  the  data  are  discussed  in  Part  2.1.  The  implied  heat 
and  momentum  sources  derived  from  Groves'  data  and  the  theoretically 
predicted  heat  and  momentum  sources  are  compared  in  Part  2.2. 

Stability  computations  for  Groves'  wind  model  are  presented  in  Part  2.5- 
In  Part  2.^  we  conclude  that  Groves'  wind  fields  are  probably  inadequate 
for  debris  cloud  advection  forecast  purposes  and  suggest  that  theoretical 
prediction  models  currently  under  development  can  be  used  to  provide 
more  reliable  results. 
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2.1  UPPER  ATMOSPHERIC  WIND  DATA 


Wind  observations  above  30  km  and  below  I50  km  are  principally 
obtained  through  rocket  based  techniques.  At  high  altitudes  a 
meteorological  rocket  releases  an  object  or  chaff  which  is  tracked 
by  radar  as  it  falls.  Lateral  motion  of  the  falling  object  then 
yields  horizontal  wind  data  and  the  local  density  of  the  atmosphere 
may  be  computed  by  observing  the  rate  of  fall.  Compared  with 
radiosonde  observations  used  below  1)0  km^  rocket  methods  are  very 
expensive  and  technologically  complex.  As  a result,  the  network  of 
rocket  launching  stations  is  quite  sparse  and  regular  observations 
are  taken  only  weekly.  Figure  2.1  shows  the  station  locations  for  the 
Meteorological  Rocket  Network  (MRN'. 

The  rocket  data  obtained  through  the  MRN  facilities  contains 
both  systematic  and  random  deviations  or  "errors".  Both  kinds  reduce 
the  usability  of  the  derived  climatological  wind  models  for  forecasting. 
We  may  further  subdivide  the  deviations  into  those  due  to  measurements 
(e.g.,  faulty  radar  techniques)  and  those  due  to  the  phenomena  (e.g., 
small  scale  eddies  inadequately  resolved  by  the  MRN  grid^.  Quiroz 
( 1969)  discusses  systematic  and  random  measurement  deviations  at 
length  and  his  findings  will  not  be  reviewed  here.  Measurement  error 
of  the  systematic  type  is  assumed  to  be  negligible,  whereas  random 
error  associated  with  the  measurements  is  assumed  to  be  removed  by 
climatological  averaging. 

Probably  the  most  obvious  source  of  systematic  deviations 
associated  with  phenomena  in  the  upper  atmosphere  is  that  produced 
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Molodeylinaza  (685,46E)  in  1969 


by  the  presence  of  tidal  winds.  Most  MRN  data  is  taken  at  local  noon. 
Thus,  if  regular  diurnal  and  semidiurnal  tidal  wind  components  are 
present,  they  will  be  interpreted  in  any  local  climatological  analysis 
as  a component  of  the  mean  wind.  Lindzen  (1967^  has  computed  the 
amplitude  of  the  tidal  winds  up  to  100  km  and  has  found  that  winds 
associated  with  the  solar  semidiurnal  and  diurnal  tides  may  be  as 
large  as  100  ms  ^ at  100  km  in  the  zonal  direction.  Measurement  of 
tidal  winds  below  60  km  shows  general  agreement  with  Lindzen 's 
computation  with  some  disagreement  evident  above  60  km  (Glass  and 
Spizzichino,  197'^^. 

Groves'  wind  model  presented  in  CIRA  (1972^  (also  Groves,  1969) 
has  been  constructed  by  grouping  MRN  and  other  data  into  monthly  or 
bimonthly  sets.  The  data  within  a set  have  been  further  subdivided 
into  four  hour  time  groups  depending  on  the  local  time  the  data  were 
taken.  The  average  within  each  group  was  computed  as  well  as  the 
mean  deviation.  Above  60  km  Southern  Hemisphere  data  were  assumed  to 
be  equal  to  Northern  Hemisphere  data.  Final  wind  model  values  were 
computed  by  an  Iterative  scheme  involving  the  average  of  the  mean 
deviations  and  the  average  of  the  group  averages  and  a weighting 
formula.  Using  an  average  of  the  group  averages  is  equivalent  in 
some  sense  to  a daily  average.  Provided  large  monthly  changes  in  the 
amplitude  and  phase  of  the  tidal  components  do  not  occur  and  data 
samples  are  present  within  each  group,  this  method  should  eliminate 
the  systematic  error  introduced  by  tides.  In  reality,  however,  many 
groups  lack  data  altogether  above  60  km  so  that  model  points  are  based 
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upon  only  one  or  two  groups  (cf.  CIRA,  1972).  We  may  conclude  then 
that  high  altitude  winds  presented  by  Groves  probably  contain  con- 
siderable bias  from  tidal  winds  superimposed  upon  the  zonal  and 
meridional  mean  winds. 

Essentially,  the  systematic  deviation  introduced  by  tidal 
components  results  from  inadequate  temporal  resolution  of  the  zonal 
mean  wind  components.  Inadequate  spatial  resolution  can  also 
introduce  systematic  deviations.  In  particular,  quasistationary 
planetary  scale  waves  as  well  as  tides  have  wind  components  which 
vary  very  slowly  over  horizontal  distances  on  the  order  of  ^,000  to 
10,000  km.  From  Figure  2.1  it  is  apparent  that  MRN  stations  are 
principally  located  in  the  northern  part  of  the  Western  Hemisphere, 
and  thus  the  network  will  be  unable  to  resolve  wind  components  associated 
with  very  long  zonal  scales. 

A comparison  of  West  European  and  North  American  data  presented  in 
CIRA  (1972^  indicates  the  presence  of  these  long  spatial  waves.  For 
example,  in  January  the  mean  zonal  wind  velocity  over  North  America  is 
~20  ms  ^ at  50  km  at  55 '^N,  while  the  mean  zonal  wind  velocity  over 
Europe  is  ms  The  difference  is  presumably  attributable  to  the 

long  wave  wind  components.  The  difference  between  North  American  and 
European  mean  zonal  wind  velocities  below  60  km  is  largest  during  winter 
and  is  consistent  with  the  observed  strength  of  planetary  scale  waves 
below  50  (van  Loon,  et  al.,  197"^).  Theoretical  calculations  of  the 
amplitude  of  stationary  planetary  waves  in  the  upper  atmosphere  indicate 
that  these  waves  may  have  sizable  amplitudes  up  to  the  mesopause  and 
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may  generate  zonal  wind  components  as  large  as  20  ms  ^ - JO  ms  ’ in 
the  mesosphere  and  lower  thermosphere  (Schoeberl,  1975)* 

Small  scale  eddies  are  probably  also  present  in  the  upper 
atmosphere  generated  by  baroclinic  instability  near  the  stratopause. 

If  their  horizontal  length  scales  are  much  smaller  than  1000  km, 
then  the  spatial  distribution  of  the  >1RN  network  will  be  inadequate  to 
properly  resolve  them.  These  eddies  will  appear  as  random  fluctuations 
in  the  MRN  data.  For  the  purpose  of  predicting  the  location  of  debris 
clouds,  these  eddies  may  be  as  important  as  the  zonal  mean  flow.  No 
information  is  available  from  Groves'  models  on  their  possible  struc- 
ture or  amplitude. 

All  of  the  phenomena  discussed  above  contribute  to  the  standard 
deviation  of  the  MRN  data  as  error.  In  Figure  2.2  the  model  values  of 
the  mean  zonal  wind  in  January  and  July  above  60  km  given  by  CIRA  (1972) 
and  the  standard  deviation  of  the  observations  from  the  model  values 
are  given.  Two  important  features  are  apparent.  First,  it  is  clear 
from  the  large  number  of  missing  standard  deviations  how  limited  the 
data  base  actually  is.  Second,  we  note  that  the  standard  deviation 
is  often  larger  than  the  mean  value  indicating  that  climatological  state 
of  the  wind  field  '^Groves'  model  ^ occurs  as  an  exception  rather  than  the 
rule. 
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2.2  HEAT  AND  MOMENTUM  SOURCES  IN  THE  UPPER  ATMOSPHERE 


The  mean  zonal  circulation  is  driven  by  external  heat  and 
momentum  sources.  In  some  cases,  these  sources  may  be  theoretically 
computed  and  a circulation  model  developed  to  compare  with  observations 
;Leovy,  196^;  Baker  and  Strobel,  1975,;  1975.  Alternatively,  a 

d D 

wind  model  derived  from  data  can  be  used  to  calculate  the  implied  heat 
and  momentum  sources  which  may  then  be  compared  to  theory  (Ebel,  197^). 
We  shall  coi.sider  the  consistency  of  computed  wind  models  and  implied 
heat  and  momentum  sources  with  their  observed  and  theoretical  counter- 
parts in  this  section. 

Leovy  (1964)  showed  that  the  westerly  stratospheric  jet  observed 
in  the  winter  hemisphere  and  the  easterly  jet  observed  in  the  summer 
hemisphere  arise  from  the  meridional  ozone  heating  gradient  in  the 
stratosphere.  Mean  zonal  wind  maximums  of  8o  ms  ^ were  computed  by 
Leovy  associated  with  mean  meridional  wind  velocities  of  0.7  ms  -. 

The  mean  zonal  wind  velocities  fluctuate  in  magnitude  throughout  the 
wintet  (Belmont,  et  al.,  1975^;  but  8o  ms  ^ is  relatively  good  agreement 
with  Groves'  (1969''  climatological  value  considering  many  of  the 
simplifications  used  by  Leovy.  However,  Groves'  mean  meridional 
velocities  are  an  order  of  magnitude  larger  than  those  computed  by 
Leovy  and  Baker  and  Strobel.  Furthermore,  their  meridional  winds  blow 
from  the  summer  pole  to  the  winter  pole,  while  Groves'  meridional  winds 
are  quite  variable  depending  upon  latitude  and  altitude. 

The  discrepancy  between  these  computations  and  Groves'  data 
may  be  due  to  several  factors.  First,  assuming  Groves  meridional  winds 
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are  correct,  the  mean  zonal  winds  (which  result  from  Coriolis  torques 
acting  upon  northward  or  southward  moving  flow)  may  be  computed  by  the 
following  equation. 

a£_s.ln.e  , 

where  Q is  the  earth's  frequency  of  rotation  and  9 is  the  latitude. 

3 is  the  Rayleigh  friction  coefficient;  v is  the  zonal ly  averaged 
meridional  velocity  of  the  wind,  and  u is  the  zonally  averaged  zonal 
velocity.  3 is  unknown  but  has  been  estimated  to  be  ~ 10  ® sec  * 

(Leovy,  196^4-^.  Using  10  ms  ^ for  v,  which  is  the  order  of  magnitude  given 
by  Groves  (I969),  gives  u = 1000  ms  which  is  inconsistent  with  the 
u values  also  given.  If  u and  v are  assumed  correct,  we  are  forced  to 
conclude  that  equation  (l)  does  not  describe  the  correct  relationship 
between  u and  v,  and  the  addition  of  a momentum  source  term,  M,  of 
unknown  value  to  the  righthand  side  of  Equation  (2.l)  is  required  to  form 
a consistent  equation  between  u and  v.  It  is  also  apparent  that  the 
magnitude  of  M must  be  quite  large.  The  presence  of  eddies  which  could 
contribute  to  M are  known  to  exist  in  winter  but  are  generally  absent 
in  summer  iKriester,  1972 However,  large  values  of  v are  also 
indicated  by  Groves  for  the  summer,  so  this  explanation  is  implausible. 

A more  complete  calculation  of  the  required  momentum  and  heat 
sources  needed  to  maintain  the  Groves  model  winds  in  the  mesosphere 
(70  - ICX)  km)  has  been  carried  out  by  Ebel  (197^)*  The  strength  of  the 
heat  sources  is  shown  in  Figure  2.5  fot  solstice  corlitions.  A 
comparison  with  the  computed  heat  sources  from  Park  and  London  ( 1975 \ 
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Fig.  2.4  — Theoretically  computed  heating  rates  for  solstice 
conditions  after  Park  and  Jordon  (1974) 


Figure  2.^,  indicates  that  the  value  of  the  heat  sources  required  to 
maintain  the  Groves  wind  field  is  roughly  an  order  of  magnitude  too 
large.  Thus,  in  agreement  with  our  above  arguments,  it  is  improbable 
that  Groves'  v values  are  consistent  with  the  u values  given. 

Equation  i>ai)  is  the  zonal  mean  momentum  equation.  If,  as 
suggested  in  Section  2.1  the  MRN  data  is  biased  by  tidal  and  stationary 
planetary  waves,  then  Equation  (l)  should  be  written  as 


§u 

^R^g  ^ ^ Sin  e = 


I a -t- 

a cos  9 &> 


(2.2) 


where  tj)  is  the  geopotential  and  the  subscript  g indicates  the  values 

given  by  Groves  (I969)  which  are  now  not  assumed  to  be  zonal  means. 

Both  tidal  and  long  wave  components  can  theoretically  produce 

meridional  wind  velocities  as  large  as  those  reported  in  Groves' 

model  (Lindzen,  19^7',  Schoebcrl,  1975^  and  in  all  probability  it  is 

these  components  that  are  reported  by  Groves.  The  larger  velocities 

permitted  by  Equation  (2.2)  arise  from  the  presence  of  a zonal  pressure 

gradient  force  on  the  righthand  side  and  the  second  term  on  the  left- 

hand  side  which  is  an  inertial  term.  Both  of  these  terms  are  much 

larger  chan  the  term  0 u for  planetary  scale  and  tidal  motions.  The 

K g 

value  of  V is  thus  not  coupled  to  u alone.  For  tidal  and  planetary 
g g 

scale  waves,  both  observations  and  theory  indicate  that  u ~v  ~10-20  ms  ^ 

in  the  stratosphere.  These  values  of  v are  more  consistent  with  the 

values  of  v and  suggest  that  the  data  are  indeed  biased  by  planetary 
g 

wave  and  tidal  components. 
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2.’. 


STABILITY  OF  MEAN  ZONAL  FLOW 


While  it  is  nearly  impossible  to  quantitatively  estimate  the 
magnitude  of  the  bias  that  long  wave  components  and  smaller  scale  eddy 
components  have  introduced  into  Groves'  wind  fields,  we  can  gain  some 
estimate  through  a stability  analysis.  Our  argument  is  as  follows: 

If  the  mean  zonal  wind  field  is  stable  to  wave  perturbations,  then 
any  finite  amplitude  eddy  disturbances  can  be  assumed  to  arise  from 
boundary  i,  tropospheric forcing.  If  the  flow  field  is  unstable,  then 
finite  amplitude  disturbances  may  arise  spontaneously  from  infinitesimal, 
local  disturbances. 

It  has  been  shown  by  Charney  and  Drazin  (1961)  that  only  large 
planetary  scale  eddies  can  propagate  into  the  upper  atmosphere. 

Synoptic  scale  disturbances  observed  in  the  troposphere  will  remain 
trapped  below  the  stratosphere.  Dickinson  (1973)  ^nd  Simmons  (1975) 
have  shown  that  the  long  wave  components  are  the  fastest  growing  modes 
for  unstable  flow  fields  similar  to  those  observed  in  the  upper  strato- 
sphere and  mesosphere.  A computation  of  the  stability  of  the  observed 
zonal  mean  flow  field  as  given  by  Groves  (CIRA,  1972)  may  thus  indicate 
where  large  amplitude  eddy  components  could  arise. 

Using  the  Charney-Stern  stability  criteria  (Charney  and  Stern, 
1962)  we  compute  numerically  the  stability  of  Groves'  mean  zonal  wind 
field.  In  the  stability  criterion  for  an  atmosphere  bounded  by  rigid 
walls,  a necessary  condition  for  instability  is  that  Q,  defined  as 


Q = 2(Q+u!') 


-2 


— ■"  5 tan  q ^ 


- sin  9 e 


32 


du) 

3z 


(2.5) 


where  z - £nvp^/p),  uu  = u/a  cos  0,  and  p is  pressure,  does  not  change 
sign  within  the  bounded  region. 

Figures  2.5  ^nd  2.6  show  the  value  of  Q computed  numerically  for 
parts  of  the  CIRA  (1972)  and  CIRA  (1965)  model  atmospheres,  respectively. 
Also  shown  is  a plot  of  the  corresponding  flow  pattern.  It  is  evident 
that  these  model  atmospheres  have  unstable  regions,  particularly  near 
the  stratopause  and  near  the  pole  at  all  levels,  as  indicated  by  the 
negative  values  of  Q.  A term  by  term  examination  of  Equation  (,2.5) 
indicates  that  the  sign  change  in  Q is  principally  a result  of  sign 
changes  in  the  last  term.  Instabilities  developing  from  this  flow 
pattern  would  thus  be  primarily  baroclinic. 

The  assumption  that  Groves'  wind  models  characterize  the  mean 
zonal  flow  field  is,  of  course,  introduced  in  this  analysis.  One 
example  where  such  an  assumption  is  certainly  incorrect  is  evident  in 
Figure  2.h  where  a patch  of  negative  Q appears  to  50 and  50  km  in  the 
wind  model  taken  from  North  American  data.  No  such  region  appears  in 
the  computations  based  upon  European  data.  The  source  of  the  sign 
change  in  Q is  the  appearance  of  a region  of  easterly  winds  in  the 
midst  of  a westerly  jet  in  the  North  American  data.  Tlie  winds  in  this 
region  are  probably  strongly  biased  by  planetary  waves  as  discussed  in 
Section  2.1  and  not  zonal  means;  hence,  the  Charney-Stern  stability 
criteria  does  not  apply. 
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llie  existence  of  high  frequency  motions  at  the  stratopause 
have  been  observed  by  Leovy  and  Akerman  1975)^  these  motions  may 

be  due  to  unstable  wind  configurations  such  as  those  given  by  the 
CIEIA  (1972,,  1965)  model  atmospheres.  In  any  event,  the  fact  that 
Groves'  models  are  unstable  indicates  that  either  eddy  components  of 
the  wind  have  biased  the  data  to  such  an  extent  that  the  wind  profile 
appears  unstable;  or,  that  eddy  components  are  present  with  sizable 
amplitudes. 
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u ms"*  0 xIO'^  sec"* 

Fig.  2.6  — Same  as  Fig.  5 for  CIRA  (1965)  model 
atmosphere  for  January  and  July 
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SUMMARY  AND  CONCLUSION 


2Jw 

The  usability  of  Groves'  climatological  wind  fields  for  fore- 
casting debris  cloud  advection  has  been  assessed  from  three  different 
viewpoints.  We  have  briefly  discussed  the  data  and  suggested  possible 
biasing  of  the  wind  field  by  tidal,  planetary  wave,  and  amall  scale 
eddy  motions.  An  examination  of  the  standard  deviation  of  the  data 
indicates  that  the  actual  structure  of  the  wind  field  in  the  upper 
atmosphere  rarely  resembles  the  climatological  data  given  by  Groves. 

We  have  also  examined  the  heat  and  momentum  sources  implied  by 
wind  fields.  The  values  of  v given  by  Groves  (I969)  are  an  order  of 
magnitude  too  large,  and  imply  heat  and  momentum  sources  much  larger 
than  expected  from  theory.  We  conclude  that  the  v values  actually 
represent  meridional  velocities  associated  with  planetary  scale  waves 
and  tides. 

Finally,  we  note  that  Groves'  zonal  wind  fields  are  unstable, 
especially  near  the  polar  stratopause.  The  instability  could  imply 
the  existence  of  large  amplitude  eddy  components  in  the  wind  field  for 
that  region.  We  conclude  then  that  Groves'  wind  fields  inadequately 
represent  the  structure  of  the  upper  atmosphere  for  the  purposes  of 
forecasting  the  advection  of  debris  clouds. 

An  alternative  to  the  use  of  Groves'  wind  fields  to  forecast 
movement  of  a debris  cloud  is  the  use  of  a theoretical  prediction 
model  which  can  be  initialized  on  a regular  schedule  or  at  the  moment 
of  debris  cloud  release.  Such  a model  is  presently  used  in  the 
lower  atmosphere  and  gives  reliable  forecasts  up  to  three  days  in 
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advance.  With  some  adaptations,  this  type  of  model  can  be  constructed 
for  the  upper  atmosphere  and  can  be  initialized  with  satellite 
radiance  data.  This  data,  which  is  in  the  form  of  temperature  fields, 
is  currently  available  up  to  60  km  (Chapman,  et  al.,  1972''  and  will 
soon  be  available  up  to  8o  km  and  higher.  Upper  atmosphere  forecast 
models  are  under  development  at  NRL  at  present.  Madala,  et  al.,  '1975) 
have  shown  that  tidal  winds  can  be  adequately  simulated  with  a spectral 
forecast  model,  and  Schoeberl  (1975)  has  been  able  to  determine  theo- 
ret  ically  the  structure  of  planetary  waves  in  the  uoper  atmosphere 
using  a similar  method. 
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Section  5 

CALCULATIONS  OF  HEATING  DUE  TO  ABSORPTION 
OF  SOLAR  RADIATION  BY  OZONE  IN  THE  STRATOSPHERE 
AND  MESOSPHERE 

L.  Baker  and  D.  Strobel 

The  thermal  structure  and  circulation  of  the  "upper  atmosphere" 
(upper  stratosphere  - mesosphere  - lower  thermosphere)  are  controlled  , 
by  the  thermodynamics  of  that  region;  specifically  the  absorption 
of  solar  radiation  and  the  radiative  transport  of  heat  energy.  As  part 
of  an  upper  atmospheric  modelling  effort  of  the  Plasma  Dynamics  Branch 
of  NRL;  a general  computer  code  was  written  to  treat  the  first  problem- 
absorption  of  insolation  by  an  absorber  with  arbitrary  distribution. 
This  code  may  be  used  through  various  entries  to  give  either  mean, 
seasonal  heatings  for  use  in  climatological  forecasts  or  to  give  point- 
by-point  local  heating  rates,  as  would  be  useful  in  tidal  calculations. 
It  can  be  useful  in  studying  the  interaction  of  circulation  and  the 
photochemistry  of  ozone  and  other  absorbers. 

Most  descriptions  of  radiative  heating  calculation  in  the 
literature  have  been  very  sketchy  as  to  the  details  of  the  integration 
procedure.  We  intend  this  report  to  fully  document  this  code.  Because 
this  code  is  not  intended  for  applications  to  the  troposphere  and  lower 
stratosphere,  we  do  not  include  scattering  or  atmospheric  refraction. 
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^.1  CALCULATION  OF  HEATING  AT  A POINT 

This  chapter  and  the  next  deal  mostly  with  geometry. 

After  determing  that  a point  is  in  fact  illuminated  by  the  sun 
(see  part  ?.2;  at  these  altitudes  scattering  is  negligible),  an  integration 
along  the  line-of-sight  is  performed  (see  2.5  for  path  length  increment 
used . 

Assume  a Cartesian  Coordinate  System  with  origin  at  the  center  of 
the  earth,  the  Z - axis  as  the  north  pole,  a vector  from  the  center 
of  the  earth  to  the  sun  lying  in  the  x-z  plane,  making  an  angel  (y 
solar  declination''  with  the  X axis  \Fig.  5*1^*  The  unit  vector  to  the 
sun  is 

S = ( cos  a,  0,  sin  O' ^ 

A point  P with  latitude  <t>  and  longitude  L is  located  by  the 
vector  P » r ( cos  cp  cos  L,  cos<J>sin  L,  sin  0 ')  where  r is  the  distance 
from  the  center  of  the  earth  and  L is  measured  from  the  X axis.  Then 
the  angle  between  the  sun  and  local  vertical  at  P,'y,  is  given  by 

P-S 

Starting  from  P we  integrate  with  step  As  along  the  line  of 
sight  to  P'  (Fig.  The  new  altitude  r'  is  calculated  by  applying 

the  law  of  cosines 

r'^  = r^  + + 2Asr  cos  y 

The  new  latitude  0'  may  be  found  from  P'  = P + AsS,  by  equating  the 
Z-components  and  solving  to  find 

0'  “ arcsin  lAs  sin  y + r sin0  ) /r' 
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Wlien  this  calculation  is  done^  using  the  colatitude  9 “ p 

(and  similarly  for  9'),  0 < S < rr,  there  is  no  ambiguity  in  the 


quadrant  of 


0'=  arc  cos  (As  sin  y + r cos  9)/r' 


Finally  we  must  find  L'.  The  easiest  way  to  do  this  (and  to  be  sure 
of  being  in  the  proper  quadrant'  is  to  use 


tan  L-  = P /P  = rw 

y X cos  4>  cos  L + As  cos  (y  /r 

sin  6 sin  L 

sin  9 cos  L + As  cos  (y)/r 

and  to  find  L'  using  the  DATAN2  function  of  the  Fortran  Library. 

Various  criteria  for  stopping  the  integration  may  be  used  (see  5-5' 
3.2  Calculation  of  Latitude,  Longitude  Limits.  Averagin)^  in 
Seasonal  Calculation. 

The  cutoff  criteria  used  below  have  multiple  purposes.  For 
the  point-wise  calculations,  they  decide  whether  an  integration 
should  be  done  or  whether  the  point  is  in  the  earth's  shadow  and 
receives  no  direct  illumination.  For  the  seasonal  calculation  they 
determine  the  duration  of  heating  (relative  to  2^  hrs.)  and  control 
the  averaging  (discussed  below^. 

A.  Latitude  Cutoff 

To  find  the  latitude  above  (or  below'  which  all  longitudes  are 


dark,  we  find  the  ray  which  grazes  the  earth  at  radius  r^  (Fig.  '^.2a'. 

For  this  ray  cos  (Fig.  5*2b,  3»2c'>;  For  illuminated  latitudes 

cos  9 is  smaller,  while  for  dark  latitudes  cos  is  larger. 
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Fig.  3.2(a)  — Basic  geometry  of  shadow  region.  Hatched  region 
is  dark.  Crosshatched  region  is  dark  all  day. 
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Fig.  3.2(b)  — Geometry  of  dark  latitudes  for  case  of 
northern  hemisphere  winter 


Fig.  3.2(c)  — Geometry  of  dark  latitudes  for  case 
of  noethern  hemisphere  summer 


'Hie  possible  magnitude  of  the  z/r  correction  to  the  cutoff 
latitude  is  not  negligible.  Consider  the  case  q/ < o (cy>  o is  completely 
analogous).  Let  A = jal.  ilien  the  cutoff  colatitude  is  given  by 

cos  (0  - a)  = ro/^r  + z) 

c o 

Then  sin(l0  - A|  ) = (1  -1/(1  + (z/r 

‘ c ' o 

For  small  z/r  , the  right-hand  side  of  this  becomes  (2z/r  )i' 
o o 

(1  + 4 — i /( 1 + z/r  or  ~ (2z/r  ^5.  For  r = 6800  km, 

* r o o o 

o 

z .X  100  km,  z/r^~.  02  sin  1 ^'.0i(  ~.2  and  | 0-A|  ^ 10°.  Note, 
however,  that  the  line-of-sight  to  points  with  colatitudes  between 
0 and  la!  lor  rr  ■ a fot  a > o''  between  the  sun  and  a point  P,  will 

pass  through  points  with  lower  altitudes.  As  the  ozone  distribution 
decreases  with  height  (except  near  mesopause,  where  the  ozone  absorption 
is  not  large),  this  effect  Is  somewhat  less  important  for  ozone  heating, 
since  such  glancing  lines  of  sight  are  usually  optically  thick  due  to 
passage  through  dense  ozone  layers  below. 

B.  Longitude  Cutoff 

Given  a partially  illuminated  latitude,  we  wish  to  inquire 
as  to  what  range  of  longitudes  are  Illuminated  (we  measure  longitude 
from  0 at  the  subsolar  point''  see  Fig.  3*5*  We  require  |R  + g 

for  all  B,  for  the  point  R to  be  Illuminated,  where  r^  is  the  radius  of 
the  earth.  Let  ~ j ^ the  radius  of  point  R.  Then  for  Q 

sin  ly  cos  0 + sin  0 cos  m cos  a,  we  can  have  the  solution  to 

|R  + gS|  =r^3=g(-  2rQ  + y''l.r''^Q'’’  - ^tA'”')  • As  B must  be  real, 

U(r''’  Q ^ - A®)>'* 
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The  limiting  tangential  case  is  for  the  equality  i.e. 

Q = + A/t.  The  negative  root  should  be  used  as  3 must  be  positive. 
Consider,  for  example,  an  equinoctal  case  a = 0,  at  the  equator  (S  = rr/2). 
Then  9 = cos  u, . = + ^/r.  Clearly  we  want  \ > 0 to  increase  the  value 

of  beyond  ^^2,  and  hence  the  negative  sign  is  mandatory. 

Averaging  and  Sensitivity  to  as. 

The  seasonal  heating  is  calculated  by  averaging  results  at  each 

latitude,  for  NAV  longitudes,  equally  spaced  from  0 v subsolar)  to  the 

maximum  Illuminated  longitude  The  daily  average  heat  input  is  then 

ffl  /tt  times  this  average.  One  may  also  use  a mean  effective  longitude, 
m 

cos  tp  - 


- %(!  cp  dcp 


sin 


in  the  calculation  of  y,  and  then  perform  one  Integration  for  each  lati- 
tude and  altitude  at  which  heating  is  desired.  Experiments  with  various 
cases  showed  that  the  results  are  not  very  sensitive  to  the  choice  of 
procedure  for  seasonal  heating.  The  differences  are  most  marked  in  the 
lower  regions  of  maximal  heating.  Here  varying  NAV  from  J to  I alters 
the  heating,  for  typical  ozone  distribution  used  \ see  below),  by  less  than 
4%.  Using  tp^  Instead  results  in  difference  of  the  same  order,  less  than 
6%;  the  larger  NAV  gave  slightly  smaller  heating  rates.  Varying  as 
from  5 to  5 4m  decreased  the  heating  rates  by  at  most  .5^,  at  the  regions 
of  maximum  heating  (less  elsewhere'). 

The  variations  above  are  for  a solsticial  case;  the  effects  of 
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altering  As  or  the  longitudinal  averaging  is  even  smaller  in  equinoctal 


cases . 

5.^  Sample  Results;  Mean  Seasonal  Heating  Due  to  Ozone  Absorption 
of  Insolation. 

The  procedures  described  above  were  used  to  find  the  solar  heating 
due  to  ozone  absorption  NAV=3^  As  = 5km  were  used, 
a^  Ozone  absorption  parameters  used. 

The  data  for  ozone  absorption  cross  section  and  solar  flux  were 
taken  from  Blake  (1973)* 

We  compared  the  calculated  values  of  "apparent  heating",  ^ 
hv  Jv  dv  , to  the  values  given  by  the  analytic  approximate  formula 
of  Lindzen  and  Will  (1973)«  Th®  agreement  is  exact  at  optical  column 
density  (cm  NTP''  u = .27«  For  larger  u our  absorption  rates  are  smaller, 
e.g.  at  u = 1 we  have  3-25  ^ 10^  ergs  * cm  ^ (cm  NTP  compared  to 
3.5  X 10^;  the  reverse  is  true  for  longer  path  lengths. 
b1  Ozone  Distribution 

For  altitudes  below  (pressures  above)  .8  mb  the  ozone  mixing 
ratios  were  taken  from  the  figures  of  Krueger  et  al  (1973^»  Table  Ic 
gives  the  pressure-latitude  distribution  of  ozone  mixing  ratio  in 
Dobson  units,  for  solstitial  and  equinoctal  northern  hemisphere  values 
are  6 mo-  different  northern  hemisphere  values.  Above  the  .8  mb 
height  level,  theoretical  ozone  distributions  were  used  v see  below'. 

1.  The  terms  "apparent  heating"  and  "actual  heating"  are  used  by 
Fukuyama  (197^b).  We  use  his  notation  in  the  integral  giving 
the  apparent-actual  heating.  Here  v is  the  frequency,  h is 
Planck's  constant,  n the  ozone  number  density,  and  J (z)  the 

photodissociation  coefficient  (l.c.,  the  product  of  solar  flux 
at  frequency  v and  height  z,  and  the  absorption  cross  section 
at  that  frequency). 


Table  la 

Ozone  Mixing  Ratio  i Dobson  Units)  as  a Function  of  Latitude  and  Pressure 

WINTER  SOLSTICE 

NORTH  LATITUDE 


0.  20.  30.  40.  60.  80.  90. 

PRESSURE 

MB. 


0.8 

2.0 

3.0 

3.5 

3.7 

6.0 

7.0 

8.0 

1.0 

k.o 

4.0 

5.0 

5.0 

6.0 

5.0 

4.0 

2.0 

9.0 

10.0 

12.0 

13.0 

13.0 

10.0 

8.0 

h.O 

lluO 

15.0 

15.0 

15.0 

16.0 

10.0 

8.0 

6.0 

16.0 

16.0 

15.0 

13.0 

12.0 

10.0 

8.0 

8.0 

16.0 

16.5 

14.0 

12.0 

11.0 

10.0 

8.0 

10.0 

15.0 

16.0 

12.0 

12.0 

12.0 

9.0 

8.0 

20.0 

11.0 

10.0 

10.0 

9.0 

11.0 

7.0 

5.0 

30.0 

8.0 

8.0 

7.5 

8.0 

10.0 

7.0 

5.0 
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Table  lb 

Ozone  Mixing  Ratio  (Dobson  Units)  as  a Function  of  Latitude  and  Pressure 

SUMMER  SOLSTICE 


NORTH  LATITUDE 


PRESSURE 

MB. 

0. 

20. 

30. 

4o. 

60. 

80. 

90. 

0.8' 

2.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

3.0 

2.0 

9.0 

9.0 

8.0 

8.0 

8.0 

8.0 

8.0 

4.0 

14.0 

14.0 

12.0 

13.0 

13.0 

8.0 

7.0 

1 6.0 

16.0 

16.0 

14.5 

13.0 

10.0 

9.0 

8.0 

1 8.0 

16.0 

16.0 

15.5 

13.0 

10.0 

8.5 

8.0 

10.0 

15.0 

16.0 

16.0 

12.0 

8.0 

8.0 

7.0 

20.0 

11.0 

10.0 

14.0 

11.0 

8.0 

7.0 

6.0 

30.0 

8.0 

9.0 

8.0 

7.0 

8.0 

6.0 

5.0 
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Table  Ic 

Ozone  Mixing  Ratio  (Dobson  Units ^ as  a Function  of  Latitude  and  Pressure 

EQUINOX 


NORTH  LATITUDE 


PRESSURE 

MB. 

0. 

20. 

30. 

4o. 

60. 

80. 

90. 

0.8 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

4.0 

1.0 

4.0 

4.0 

5.0 

4.0 

6.0 

7.0 

8.0 

2.0 

9.0 

9.5 

8.0 

10.0 

11.0 

7.0 

6.0 

4.0 

13.5 

14.0 

12.0 

13.0 

10.0 

6.0 

4 .0 

' 6.0 

16.0 

16.5 

14.0 

12.0 

8.0 

5.0 

3.0 

8.0 

16.5 

16.5 

14.5 

11.0 

8.0 

5.0 

2.0 

j 10.0 

16.0 

16.5 

14.5 

10.0 

8.0 

4.0 

1.0 

20.0 

11.0 

12.0 

12.0 

8.0 

8.0 

3.0 

1.0 

1 30.0 

8.0 

3.0 

7.7 

8.0 

8.5 

8.2 

7.5 

Heating  functions  were  computed  with  the  infrared  radiative  cooling 


parameterization  of  Dickinson  (1975)«  This  approximation  becomes  poor 


above  a 65  km  but  nothing  better  is  currently  available.  Below  50  km 


the  relaxation  rate  is  taken  as  fixed  and  equal  to  the  rates  at  50  km. 


The  cooling  1,  to  space 'i  term  is  taken  to  behave  as  c(70)  exp  [-(Z-70)/3] 


z in  km. , above  70  km. 


The  relaxation  rate  due  to  ozone  photochemistry  is  an  additional 


contribution  to  the  thermal  relaxation  rate.  It  is  calculated  following 


Leovy  (1964).  Comparison  with  the  work  of  Blake  1, 1970)j  however  a more 


detailed  reaction  rate  set,  is  quite  favorable.  This  contribution  to  k , 

R 

K , is  proportional  to  the  local  ozone  heating.  Then  k = K + K , 

'T  R 03  IR 

K found  from  Dickinson  (1975^  described  above. 

I R 

The  ozone  heating  is  found  using  the  program  discussed  above. 


Below  50  km  the  ozone  distribution  used  is  that  of  Kreuger  et  al.  1 1972 ^ . 
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above  50  km,  the  values  are  from  Crutzen  1771)  for  the  summer  and 
Fukuyama  (197^a'  for  the  winter.  For  equinox  a latitudinally  independent 
distribution  that  was  intermediate  between  winter  and  summer  values  was 
used.  Contours  of  these  distributions  are  shown  in  Figs.  5.2a  and  5.5a. 

In  computing  the  heating  rate  above  70  km  the  heating  function 

Fukuyama  197^t>)  is  used  and  is  assumed  proportional  to  n ■ hv-E-D) Jvdv 

where  E,D  are  excitation  and  dissociation  energies.  Below  70  km  we  use 

an  "apparent"  heating  rate,  assuming  that  the  chemical  energy  is  locally 

converted  to  thermal  energy.  Clearly  the  transition  between  these 

limiting  cases  is  not  so  abrupt,  but  this  is  adequate  for  our  calculation. 

The  resultant  heating  rates  are  shown  in  Figs  5-5h.  Note  the  "island" 

of  heating  in  the  upper  low- latitude  winter  hemisphere  due  to  large 

Og  density.  There  is  a relative  maximum  of  heating  at  the  summer  pole 

~35  km  which  is  not  seen  due  to  the  contour  interval.  Combined  with 

an  assumed  T (z)  (COSPAR  196l^  the  infrared  transfer  parameterization 
o 

described  above  gives  the  net  heating/cooling  rates  shown  in  Figs  5.^c, 

3.5c.  The  thermal  relaxation  rates.  Figs.  5.^d,  5.5d  are  also  found 
as  described  above.  Using  k^  and  the  net  heating  Q,  a radiative  equili- 
brium temperature  may  be  found,  as  in  Figs.  3.‘'e,  Ihe  solstitial 

results  are  closer  to  the  profile  assumed  by  Trenberth  1975'',  than  that 
of  Leovy  (196^a'i;  Leovy's  results  evidence  much  colder  winter  pole  values. 
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Fig.  3.4(a)  — Contours  of  ozone  number  density  assumed-solstice 
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P’ig.  3.4(d)  — Thermal  relaxation  rate,  day 
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Kin.  3.4(e)  — Radiative  equililmum  temperatures 
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Fig.  3.5(a)  — Contours  of  ozon**  numlwr  density  assumed  - equinox 
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Fi«.  3.5(b)  — ResulUnt  hoatinR  rate  due  to  insulation  in  dep-ees  day 
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F\^.  3.5(e)  — Radiative  equilibrium  temperatures 


ion  of  the  Heating  Function  & Models  of  the  Mesospheric 
Circulat ion. 

In  a previous  report  Baker  and  Strobe  1,  1975;  hereafter  referred 
to  as  Report  we  constructed  linear  separable  models  to  describe  the 
mean  circulation  of  the  upper  stratosphere  and  mesosphere.  This  circu- 
lation can  have  a profound  effect  on  radioactive  debris  patches  created 
by  high  altitude  nuclear  explosions  ' Zalesak  and  Coffey,  1975)*  The 
subsequent  beta  decay  results  in  long  lasting,  widespread  ionization 
that  can  severely  affect  the  performance  of  communication  systems.  Ac- 
curate prediction  of  system  performance  depends  critically  on  our  abil- 
ity to  describe  the  movement  of  debris  patches  by  mesospheric  wind 
systems,  of  which  the  mean  circulation  is  one  important  component. 

Accurate  simulation  of  the  mean  circulation  in  the  upper  strato- 
sphere and  mesosphere  which  is  primarily  driven  by  differential  radiative 
heating  due  to  ozone  absorption  requires  accurate  parameterization  of 
radiative  heating.  In  Report  I the  emphasis  was  on  the  importance  of 
other  parameters,  particularly  the  Newtonian  (radiative'  cooling  rate 
and  the  Rayleigh  friction  coefficient.  The  radiative  heating  functions 
of  Leovy  f l^Gha^  were  used  in  the  calculation. 

In  this  report  we  describe  calculations  with  the  improved  heating 
functions  described  above  and  the  Dickinson  (1975^  parameterization  of 
infrared  radiative  cooling.  Our  improved  models  accurately  simulate 
the  observed  mean  zonal  mesospheric  circulation  for  both  solst’tial 
and  equinoctial  conditions. 

In  the  course  of  our  research  we  have  discovered  a number  of 
inconsistencies  in  the  Leovy  (T)^^-a.b)  heat  functions.  Briefly  the 
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inconsistencies  of  Leovy  il96^<a,b'  relate  to  the  assumed  equilibrium^ 
basic  state  of  the  atmosphere  (mean  density,  pressure,  temperature,  and 
wind  fields^  and  the  corresponding  differential  heating  functions  which 
are  temperature  dependent.  For  example,  any  latitudinal  variation  in 
the  equilibrium  basic  temperature  field  requires  a corresponding  zonal 
geostrophic  wind  to  describe  the  equilibrium  state  of  the  atmosphere. 

This  zonal  geostrophic  wind  was  ignored  in  some  of  Leovy's  models. 

The  ozone  heating  rate,  Qvrr,cp);  where  tt  is  the  vertical  coordinate 
log  pressure  coordinates,  see  Report  l'  and  ^ is  latitude,  is  computed 
for  an  assumed  ozone  density  distribution  as  described  above.  The  in- 
frared radiative  cooling  rate  is  based  on  Dickinson's  1973)  results. 

The  thermal  relaxation  rate,  kj^(rr,cp^  is  given  by  the  sum  of  the  IR 
cooling  rate  and  the  ozone  photochemical  damping  rate  calculated  by 
Leovy  (196^J-b)  and  is  consistent  with  the  results  of  Blake  (1973)* 

Given  the  functions  Q(rr,«p^  and  (rr^tp)  we  can  compute  the  lati- 

tudinally  averaged  values  of  k , denoted  by  <k  \tt^>  and  expand  the 

R R 

heating  function  in  a modal  representation 


Q(rr,ep) 


(3.1) 


m 

as  defined  in  Report  I and  where  y - sin  cp.  We  calculate  Sm.TT^  by  the 

usual  technique  of  assuming  the  expansion  form  (l'  and  use  the  ortho- 

dT] 

rn 

gonality  of  the  — functions.  These  functions  were  evaluated  at  each 
TT  level  of  the  model.  'Pwo  modes  were  required  for  solstitial  condi- 
tions whereas  only  one  mode  sufficed  for  equinoctial  conditions.  Their 
numerical  values  are  shown  in  Figs.  '^.6  and  7.7  respectively. 


■ * . ■>.  I"  " 


For  equinoctial  conditions  the  latitudinal  average  of  k (tr,^) 

R 

weighted  by  cos  co  to  account  for  surface  area  variation  with  latitude 
was  used  for  <kj^  (rr'>.  For  solstitial  conditions  it  is  necessary  to 
insure  the  correct  total  thermal  forcing  that  drives  the  circulation. 

In  radiative  equilibrium 

Teq  = To(n)  + 

R 

where  Tgq  is  the  radiative  equilibrium  temperature  and  Tq  its  mean 
value.  The  pole  to  pole  temperature  difference 

ATeq  = rQ(TT,T90°)  - Q(TT,-90°)]/<kj^(TTi> 
determines  <k^(TT^>.  These  values  for  <kj^(rx)>  were  slightly  smaller 
(<1G^)  than  the  "weighted"  mean  kj^(n^fp)  cos  tp.  In  addition  the  weigiited 
arithmetic  <kj^  cos  > and  harmonic  <(kj^  ^cos  cp  '>  ^ means  differ  only 
slightly. 

The  mean  temperature  profile  To(tt'  and  the  corresponding  static 

dTo 

stability  parameter  R (tts  ^ function  of  To(n^  and  — r-  , must  be  consis- 

S (ITT 

tently  selected.  Based  on  the  CIRA  ^1961^  mean  temperature  profile^  the 
corresponding  height  profile  of  Rg'rr^  is  illustrated  in  Fig.  '^8.  Ihe 
Rayleigh  friction  coefficient,  R , was  taken  to  be  1.5  x 10  ” sec  ^ in 

r 

accord  with  other  models. 

The  model  results  for  solstitial  conditions  are  shown  in  Figures 
3. 9a  -3.9d.  We  note  from  Fig.  3*9®  that  the  model  predicts  zonal  jets  of 
approximately  correct  velocities  and  at  the  observed  tieights.  For  ex- 
ample Murgatroyd  (1970^  gives  observed  peak  velocities  of  m sec 
for  the  winter  jet  and  -60  m sec  for  the  summer  jet.  'Ihe  winter  jet 
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is  centered  approximately  10  km  below  the  summer  jet  and  they  are  not 
on  the  same  geopotential  surface.  We  expect  this  result  due  to  the  ad- 
ded symmetric  heating  mode  which  forces  an  asymmetric  zonal  velocity 
field.  The  fact  that  the  resultant  zonal  velocity  field  is  in  good 
agreement  with  observations  would  suggest  that  a symmetric  component  of 
heating  is  an  essential  element  of  mesospheric  circulation  at  solstice. 

At  solstice  the  winter  mesosphere  is  predicted  to  be  approximately 
isothermal  in  agreement  with  observations  ,Fig.  No  attempt  was 

made  to  accurately  model  the  lower  stratosphere,  which  is  controlled  by 
tropospheric  forcing.  Thus  the  predicted  temperature  minimum  at  the 
winter  pole  is  Ji0‘^  K colder  than  observed,  whereas  at  the  summer  pole 
the  difference  is  only  5°  K (Murgatroyd,  1970). 

The  meridional  and  vertical  velocity  fields  are  small  cf.  Figures 
3.9c  and  3.9<^^*  As  discussed  in  Report  I they  are  very  sensitive  to  ihe 
magnitude  of  the  Rayleigh  friction  coefficient.  Tlie  vertical  velocities 
(Fig.  ■^'.9d''  show  sinking  motion  at  the  winter  pole  at  twice  the  rate  of 
the  rising  motion  at  the  summer  pole  due  to  the  symimetric  heating  mode. 
Also  the  meridional  velocity  field  exhibits  similar  asymmetries  due  to 
this  heating  mode.  We  anticipate  that  the  inclusion  of  non-linear  ad- 
vection  would  result  in  reduced  meridional  and  vertical  velocities  at 
the  winter  pole. 

At  equinox  the  model  results  are  symmetrical  with  respect  to  the 
equator  and  only  one  hemisphere  is  illustrated  in  Figs.  . 10a~3 . lOd . The 
zonal  jet  velocity  is  comparable  in  magnitude  to  solstitial  conditions, 
but  are  accompanied  with  peak  meridional  winds  only  one-half  as  strong 
as  at  solstice.  The  poleward,  symmetric  about  the  equator,  motion  re- 
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suits  in  a "2  - cell"  circulation  rather  than  the  predominate  1 - cell 
circulation  at  solstice.  At  midlatitudes  the  equinoctial  results  for 
zonal  wind  and  temperature  are  in  good  agreement  with  observations. 
However  at  the  equator,,  where  the  quasi-geostrophic  assumption  breaks 
down,  there  are  observed  zonal  jets  ~ f|-0  m sec  ^ not  predicted  by  our 
mode  1 . 

5 . 5 Conclusion 

Linear  zonally  symmetric  models  can  accurately  simulate  the  ob- 
served zonal  wind  and  temperature  fields  particularly  at  mid  and  high 
latitudes.  Based  on  realistic  values  for  the  Rayleigh  friction  co- 
efficients we  predict  peak  vertical  velocities  of  0.^  - 0.6cm  sec  - at 
the  poles  during  solstice  and  equinox.  The  peak  mean  meridional 
velocities  occur  at  midlatitudes  and  are  < Im  sec  . 

It  should  be  noted  that  only  ozone  heating  was  included  in  our 
model.  The  next  step  is  to  investigate  the  importance  of  atmospheric 
heating  due  to  molecular  oxygen  absorption  in  the  Schumann-Runge  bands 
above  ~ 70km  on  mesospheric  circulation.  In  addition  the  non-l.TE  ef- 
fects on  the  infrared  radiative  cooling  rate,  neglected  in  the  Dickinson 
(1977)  parameterization  must  also  be  evaluated.  These  sources  and  sinks 
of  heat  energy  tend  to  offset  each  other.  We  do  not  anticipate  that  the 
inclusion  of  tliese  processes  will  dramatically  alter  the  results  pre- 
sented in  this  report.  At  equinox  the  physics  of  the  equatorial  jets 
must  be  included  for  accurate  simulation  in  the  equatorial  regions  of 
the  mcsosphi're. 
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Section  I 

A SEMI-SPECTRAL  NUMERICAL  MODEL  FOR  FORCED, 

X'ERTICALLY  PROPAGATING  PLANETARY  WAVES 

Part  I-Application  of  the  Model  to  Linear  Diurnal 
and  Semi-Diurnal  Atmospheric  Thermal  Tides 

R.V.  Madala,  S.A.  Piacsek,  and  S.T.  Zalesak 
The  study  of  non-linear  interactions  of  forced  planetary  waves 
with  one  another  and  with  the  mean  current  can  best  be  treated  as  a time 
dependent  initial-value  problem.  However,  this  approach  suffers  from  the 
disadvantage  that  some  of  the  free  modes  of  the  atmosphere  will  be  excited 
besides  the  desired  response  to  a particular  forcing,  due  to  inaccurate 
initial  conditions.  As  shown  by  Lindzen  et  <"1.  (1968)  the  numerical 
models  also  suffer  from  another  disadvantage  in  that  the  vertically 
propagating  modes  are  reflected  back  by  a rigid  or  inaccurate  upper 
boundary,  leading  to  change  in  the  structure  of  these  waves. 

Yanowitch  (I969)  and  Houghton  and  Jones  I96 ' have  investigated 
the  vertical  propagation  of  gravity  waves  in  isothermal  atmospheres,  and 
studied  the  reflection  and  absorption  of  these  waves  by  an  exponentially 
increasing  friction  profile  and  Rayleigh  viscosity,  respectively.  Hong 
and  Lindzen  1973)  have  simulated  the  tides  numerically  in  the  presence 
of  a geostrophic  mean  wind,  but  obtained  convergent  results  only  for 
the  semidiurnal  tide.  In  each  of  these  studies,  only  sinusoidal  waves 
of  the  form  were  considered,  leading  to  an  eigenvalue  problem  of 

ordinary  differential  equations  in  the  height  z.  The  gravity  wave  studies 
were  concerned  only  with  plane  geometry,  whereas  the  tidal  studies  in- 
volved spherical  harmonics  and  treatment  of  the  poles;  on  the  other  hand, 
the  latter  did  not  include  the  formation  of  critical  layers  i.c.,  where 
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the  phase  speed  of  the  wave  equals  the  wind  speed  at  that  level),  because 
both  the  2h  and  12  hour  tides  propagate  around  the  Earth  much  faster 
than  any  wind  can  move. 

The  main  thrust  of  the  present  paper  is  to  utilize  the  informa- 
tion derived  by  the  foregoing  studies  on  vertical  propagation  and 
reflection,  and  to  consider  the  alternate  method  of  asymptotically 
reaching  a quasi-steady  oscillatory  state  of  forced  motions  via  initial 
value,  time-dependent  methods.  This  seemed  to  be  desirable  not  only 
to  avoid  the  difficulties  associated  with  inverting  the  diurnal  matrix, 
but  to  formulate  more  general  techniques  for  non-linear  problems.  The 
tides  were  deemed  an  ideal  problem  to  test  th.e  linear  version  of  the 
model,  since  both  analytical  and  numerical  at  least  for  semi-diurnal 
tide*  solutions  are  available. 

The  disadvantages  mentioned  earlier  are  circumvented  in  this 
model  by  introducing  a transient  Rayleigh  friction  decaying  in  time) 
to  get  rid  of  the  unwanted  free  modes,  and  a sponge  layer  below  the 
upper  boundary  to  absorb  the  vertically  propagating  modes.  As  a test 
case,  the  linear  version  of  the  model  is  applied  to  the  semi-diurnal 
and  diurnal  thermal  tides  in  the  lower  and  middle  atmospheres  (<100  km'. 
The  numerical  results  are  compared  with  the  analytical  results  obtained 
by  Lindzen  (1971'* 

Since  ve  are  dealing  with  a small  number  of  planetary  waves, 
a semi-spectral  model  with  Fourier  expansions  in  the  E-W  direction  is 
more  efficient  than  a ~^-D  grid  point  model  because  only  a limited 
number  of  modes  arc  needed  to  describe  these  waves.  In  general,  at 
least  9-10  mesh  points  per  wave  length  are  needed  to  describe  a 


sinusoidal  wave  accurately,  so  for  phenomena  containing  a narrow  spectral 
baud  the  saving  in  computer  storage  can  be  very  large.  }Iowever, 
spectral  models  are  difficult  to  use  in  a system  with  height  as  the 
vertical  coordinate,  since  the  corresponding  pressure  tendency  equation 
contains  dependent  variables  in  the  denomenator.  Therefore,  in  the 
present  model,  a quantity  tt  = - H^Xn  P/P^  is  used  as  the  vertical  coordin- 
ate instead  of  height.  The  spectral  equations  in  the  rr-coordinate 
system  are  less  complicated  than  those  in  the  z-coordinate  system. 

The  physical  model  (described  in  Section  'hi)  used  for  the  test 
case  is  the  same  as  the  one  described  by  Lindzen  (1971).,  except  for  n 
being  the  vertical  coordinate  instead  of  height.  The  numerical  techniques 
are  described  in  Section  4.9.  Ilie  model  results  are  compared  with 
Lindzen 's  (1971)  results  in  Section 
1 . 1 'fhe  Physical  Model 

Following  previous  modelers,  we  shall  make  the  following 
assumptions ; 

a.  The  atmosphere  is  a perfect  gas  in  local  thermodynamic 
equilibrium,  with  the  gas  constant  uniiorm  in  height. 

b.  ITie  earth  is  a perfect  sphere  with  no  topography, 

c.  The  atmosphere  is  thin  compared  to  the  earth's  radius. 

d.  Gravitational  potential  due  to  moon  can  bo  neglected, 

e.  The  basic  (or  background)  atmosphere  is  in  hydrostatic 
balance  with  no  motion, 

f.  Tidal  oscillations  are  in  quasi-hydrostatic  balance,  and 

g.  Equations  of  motion  for  the  tidal  oscillations  can  be 
linearized. 


Lindzen  1971)  gave  some  justification  of  making  these  assumptions 
in  a study  on  atmospheric  thermal  tides;  for  more  details,  the  reader 
is  referred  to  this  monograph. 

With  these  assumptions,  the  relevant  equations  (see  Appendix  A 

for  list  of  symbols)  governing  the  tidal  oscillations,  in  a system  with 

TT  = - H P/P  as 
o o 
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Equations  (4.1'  and  (4.2)  represent  momentum  conservation  in 

the  E-W  and  N-S  directions  respectively.  The  terms  F and  F are  frictional 

u V 

terms  and  are  intended  mainly  to;  (a)  get  rid  of  the  unwanted  free 
inertia-gravity  modes  generated  by  the  imperfect  initial  conditions;  and 
b)  avoid  reflections  of  the  vertically  propagating  waves  from  the  top 
boundary.  Equations  1,4.5)  and  (4.4)  represent  the  hydrostatic  balance  and 
mass  continuity,  respectively,  and  (4.5)  is  the  thermodynamic  equation. 

The  term  Q on  the  right  side  of  (4.5)  represents  diabatic  heating  and 
contains,  in  this  model,  absorption  of  radiation  by  water  vapor  and 
ozone.  The  values  of  Q are  computed  by  Leovy  (1964)  for  various  seasons 
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of  the  year,  and  their  diurnal  and  semi-diurnal  components  in  terms  of 
Hough  modes  are  given  by  Lindzen  (1971)*  The  latter  values  of  Q are 
used  in  this  model.  Equation  (h.6)  is  the  tendency  equation  for  the 
geopotential  of  the  lower  boundary. 
h.2  The  Numerical  Model 

Unlike  in  the  usual  spherical  harmonic  (spectral'  models^  in 
the  semi-spectral  model  the  variables  are  expanded  in  Fourier  series 
in  the  E-W  direction  only,  and  are  specified  at  discrete  grid  points 
in  the  N'-S  and  vertical  directions.  The  equations  governing  the 
coefficients  of  the  Fourier  modes  are  obtained  by  substituting  the 
Fourier  expansions  of  all  the  dependent  variables  in  Equations  (4.1' 
through  (4.6).  The  resultant  partial  differential  equations  in  latitude, 
height,  and  time  are  solved  by  finite  difference  methods. 

The  distance  between  the  two  poles  is  divided  into  equally 
spaced  mesh  intervals  with  a grid  spacing  of  2.5^^  latitude.  The  thermo- 
dynamic variables  and  vertical  velocity  are  specified  at  one  set  of 
grid  points  (Figure  4.1''  and  the  horizontal  components  of  the  velocity 
are  specified  at  another  set  located  half-way  between  them.  The 
vertical  direction  is  divided  into  a number  of  layers  (Figure  4,2^. 

All  the  variables  except  the  vertical  velocity  are  specified  at  the 
middle  of  each  layer,  and  the  vertical  velocity  is  specified  at  the 
boundaries  separating  the  layers.  Below  100  km  a layer  thickness  of 
1 .on  is  used  for  the  diurnal  and  2 km  for  semi=diurnal  tide,  since  the 
former  exhibits  a more  rapid  vertical  variation.  Above  100  km  the 
respective  thicknesses  are  2 and  4 km,  since  the  Rayleigh  friction 
damping  increases  rapidly  with  height. 
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Fig.  4.2  — Grid  system  in  the  N-S  direction 
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Time  Integration  Scheme 


The  model  is  integrated  by  using  a two  step  predictor-corrector 
(leapfrog-trapezoidal'  scheme.  When  applied  to  a wave  equation, 
this  scheme  has  a neutral  amplitude  growth  with  a very  small  phase 
error  (Kurihara,  1965^.  The  two  steps  of  integration  for  the 
equation  ^ = F H'  are: 

11*  = H*"  + 2Dt  F H*")  Leapfrog  Predictor  (^.7) 

^ |F'H  ) + F(H^''|  Trapazoidal  Corrector  (^.8) 

where  the  superscript  represents  the  time  level  and  H is  the  value  of 
H predicted  by  leapfrog  scheme.  Even  though  the  scheme  generates  a 
computational  mode,  the  latter  never  grows  big  enough  to  cause  appreciable 
time  splitting. 

The  Sponge  Layer  and  Rayleigh  Friction 

One  of  the  main  disadvantages  of  using  the  initial-value  problem 
approach  for  the  tidal  problem  is  that  many,  if  not  all,  of  the  free 
modes  of  atmospheric  oscillations  will  be  excited  due  to  imperfectly 
posed  initial  conditions.  These  are  eliminated  in  this  model  by  intro- 
ducing a transient  Rayleigh  friction  of  about  I.5  x 10  sec  ^ (which 
corresponds  to  an  e-folding  time  of  about  six  days  for  atmospheric 
oscillations'  throughout  the  lower  100  km  of  the  atmosphere.  It  is 
gradually  reduced  over  a period  of  25  days  to  a value  of  about  10  ^sec  ^ 
(yielding  an  atmospheric  e-folding  time  of  I08  days ^ . Since  the  final 
value  of  the  Rayleigh  friction  is  neglegibly  small,  the  results  obtained 
from  the  model  will  approximately  correspond  to  those  of  a non-viscous 
atmosphere. 
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For  an  atmosphere  with  little  or  no  friction,  the  initial  value 
approach  suffers  from  another  disadvantage,  namely,  that  the  vertically 


propagating  waves  are  reflected  back  from  the  top  boundary  and  contaminate 
the  solution.  Tliese  reflections  are  eliminated  in  the  model  by  introduc- 
ing a sponge  layer  with  a thickness  between  50  km  to  80  km  above  the 
region  of  interest.  In  this  layer,  the  Rayleigh  friction  is  increasing 
with  height,  thereby  absorbing  the  vertically  propagating  oscillations 
as  they  move  up  (and  down)  through  the  region.  This  upper  layer  friction 
remains  constant  in  time;  in  order  to  maintain  continuity  with  the 
diminishing  lower  layer  friction,  we  must  lower  the  boundary  of  the  lower 
and  upper  regions  appropriately  with  time,  approaching  80  km  near  the 
end.  The  terms  and  F^  in  Equations  (^.1  and  (h.2''  take  the  following 
form: 

F = - K (tt)  u (I4.9) 

and  F = - K ; tt  ’ V ( 4 . 10 ' 

V z 

where  (rr''  =0.25  x 10  exp  (n  + Ttp  - 95  V?  tTt 

=K  (tT-i  ) TT„  <TT<TT,  ] 

Z 1 B-  - j 

Figure  illustrates  the  behavior  of  K vTT^rr^(t)).  ! 
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i 

Boundary  and  Initial  Conditions  I 

i 

In  order  to  solve  Equations  (It.l)  through  (ii.6''  numerically,  it  is 
sufficient  to  specify  the  values  of  the  geopotential  and  the  vertical 
velocity  at  the  two  poles,  and  vertical  velocity  at  top  of  the  model. 

All  of  these  are  assumed  to  be  zero.  The  absence  of  real  friction 
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and  advection  makes  these  simple  conditions  possible. 


Initially^  all  the  tidal  perturbations  are  assumed  to  be  of 
zero  amplitude.  The  basic  atmosphere  is  in  hydrostatic  balance  with 
no  motion.  The  basic  atmosphere  for  the  diurnal  tide  is  assumed  to  be 
isothermal  with  a temperature  of  26o°K.  The  temperature  distribution 
for  the  basic  atmosphere  in  semi-diurnal  case  is  given  in  Figure 
In  the  lowest  100  km^  this  is  the  same  as  the  equatorial  temperature 
profile  given  by  Lindzen  (1971^« 

The  diabatic  heating  is  that  of  Lindzen  (1971  and  consists  of 
linear  combination  of  the  (2,2),  {2,h),  and  \2,6)  Hough  modes  for  the 
semi-diurnal  tide^,  and  the  '1,-4),  1, 1.-2),  1,1),  (l,3)j  and  (1,5)  Hough 

modes  for  the  diurnal  tide.  For  the  semi-diurnal  tide,  all  the  three 
Hough  modes  propagate  vertically  for  the  temperature  distribution  used 
in  the  calculations  (Fig.  4.3^  except  the  (2,2)  mode  between  50  km  and 
3o  km  altitude.  For  an  isothermal  basic  temperature,  the  modes  (1,-4'' 
and  ;l,-2'  are  trapped  completely,  while  the  other  three  modes  propagate 
vertically  everywhere.  The  trapped  modes  dominate  at  higher  latitudes 
while  the  propagating  modes  dominate  at  latitudes  between  the  equator 
and  30^^.  The  modes  (2,6)  and  (1,5)  have  the  shortest  vertical  wave 
lengths  for  semi-diurnal  and  diurnal  tides,  respectively,  with  wave 
lengths  of  about  33  hm  and  7 km.  The  reader  is  referred  to  Lindzen  (1)71 
for  more  details  about  the  structure  and  behavior  of  these  modes. 

The  Diurnal  Tide 

In  order  to  see  how  well  the  model  reproduced  both  the  trapped 
and  propagating  modes,  a comparison  is  made  with  Lindzen's  (1)71)  results 
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Fig.  4.3  — Mean  atmospheric  temperature  (°K) 
for  the  semi-diurnal  tide 
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at  75°j  ^5°  and  15°  latitudes.  Amplitudes  and  phases  of  the  three 
components  of  the  tidal  velocity  are  given  in  Figures  ^.5  through  ^.10  at 
75°  latitude,  in  Figures  ^4^.  11  through  k.lG  at  ^5°  latitude,  and  in  Figures 
4.17  through  4.22  at  15°  latitude.  It  can  be  seen  from  these  figures  that 
the  model  reproduced  the  analytical  results  at  75°  and  45°  latitudes 
very  well.  However,  the  numerical  results  deviated  slightly  from  the 
analytical  results  at  15°  latitude.  This  can  be  attributed  to  poor 
vertical  resolution,  insufficient  to  resolve  the  ;1,5)  Hough  mode. 

According  to  Gramme Itvedt  (I969),  in  general  one  needs  more  than  I5 
grid  points  per  wave  length  to  reproduce  the  amplitude  and  phase  of  a 
propagating  wave.  With  1 km  grid  length,  we  have  only  about  7 grid 
points  to  resolve  the  (1,5^  Hough  mode,  therefore,  one  can  attribute 
the  observed  deviations  to  truncation  errors.  It  is  important  to  note 
that  the  vertical  velocities  presented  in  Figures  4.9,  4.15  and  4.21  repre- 
sent rather  than  In  order  to  achieve  a comparison,  the 

df^  dt  dt 

values  corresponding  to  Lindzen's  z calculations  have  been  computed. 
Semi-Diurnal  Tide 

It  is  easier  to  simulate  the  semi-diurnal  tide  numerically  than 
the  diurnal  tide,  since  the  fo’-mer  has  a less  complicated  vertical 
structure  than  the  latter.  To  show  the  accuracy  of  the  model  in  simu- 
lating the  semi-diurnal  tide,  the  amplitudes  of  u at  75  latitude,  v 
at  45°  latitude,  and  w at  15°  latitude  are  plotted  against  the  analytical 
results  in  Figures  4.23,  4.24,  and  4.25  respectively.  It  can  be  clearly 
seen  from  the  figures  that  the  model  reproduces  the  analytical  semi-diurnal 
tide  almost  exactly. 
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4.6  — Phase  of  the  westerly  velocity  of  the 
diurnal  tide  at  75°  latitude 
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Fig.  4.7  — Amplitude  of  the  northerly  velocity  (cm  see""'  ) 
of  the  diurnal  tide  at  75°  latitude 
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Fig.  4.8  — Phase  of  the  northerly  velocity  of  the 
diurnal  tide  at  75°  latitude 
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Fig.  4.9  — Amplitude  of  the  vertical  velocity,  dw/dt,  of  the 
diurnal  tide  at  75°  latitude  (units;  cm  sec”'  ) 
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Fig.  4.10  — PhasE'  of  the  vertical  velocity  of  the 
diurnal  tide  at  75°  latitude 


10  10*  I o'  I o’ 


amplitude  - ZONAL  VELDCllr  ICM/SECI 

Fig.  4.11  — Amplitude  of  the  westerly  velocity  (cm  sec~* ) 
of  *he  diurnal  tide  at  45°  latitude 
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Fig.  4.12  — Phase  of  the  westerly  velocity  of  the 
diurnal  tide  at  45°  latitude 
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Fig.  4.14  — Phase  of  the  northerly  velocity  of  the 
diurnal  tide  at  45°  latitude 
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Fig.  4.17  — Amplitude  of  the  westerly  velocity  (cm  sec 
of  the  diurnal  tide  at  15°  latitude 
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Fig.  4.18  — Amplitude  of  the  northerly  velocity  (cm  sec“*) 
of  the  diurnal  tide  at  15°  latitude 
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Fig.  4.23  — Amplitude  of  the  westerly  velocity  (cm  sec~^) 
of  the  semi-diurnal  tide  at  75°  latitude 
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Fig.  4.25  — Amplitude  of  the  vertical  velocity,  dw/clt,  of  the 
semi-diurnal  tide  at  15°  latitude  (units;  cm  sec~l) 
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Section  'j 


BETA  PATCH  DEVELOPMENT 
S.Zalesak;  J.  Blocks  D.  Strickland 
ADVECT  is  a Lagrangian  computer  program  which  tracks  the  advection 
of  a passive  nuclear  debris  cloud  by  following  the  motion  of  nets  of 
tracer  particles  defining  the  exterior  and  interior  surfaces  of  the  cloud. 
The  wind  field  responsible  for  the  advection  is  input  as  an  external 
parameter.  For  the  present  study  the  wind  model  used  is  the  empirically 
based  one  given  by  Groves  (I969).  CIRA  (1972)  in  the  form  of  mean  seasonal 
longitudinally  averaged  zonal  and  meridional  fields  for  latitudes  from 
80°  S to  80*^  N and  altitudes  from  25  km  to  IJO  km. 

When  ADVECT  was  applied  to  typical  nuclear  debris  clouds  for 
times  up  to  ^8  hours  after  release , results  differed  significantly  from 
those  of  other  models,  specifically,  the  one  used  by  the  WEPH  code. 

WEPH  produces  a cloud  expanding  uniformly  in  the  horizontal  direction 
and  at  all  times  centered  at  the  release  point:  ADVECT  predicts  elongated, 

almost  thread  like  clouds,  with  a predominately  east-west  orientation 
with  generally  significant  displacements  of  the  center  of  mass  from  the 
release  point.  Since  beta  patch  patterns  are  similar  in  shape  to  the 
debris  clouds ^it  is  likely  that  the  differences  in  cloud  advection 
between  the  two  models  would  also  be  reflected  in  the  predictions  of 
communications  interference. 

'fhe  communications  environment  at  any  given  time  after  a nuclear 
event  is  characterized,  as  far  as  beta  effects  are  involved,  by  the 
distribution  of  free  electrons  in  the  atmosphere.  To  arrive  at  this 
distribution  it  was  necessary  to  add  to  ADVECT  subroutines  which  would 
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generate  a properly  distributed  beta  flux,  calculate  the  resulting 


i 

I 

ionization  rates  for  regions  conjugate  to  the  clouc^  and  finally  produce  j 

a quasi-equilibrium  density  distribution  for  free  electrons.  The  analytic 
model  of  Knapp  and  Fischer  (1970)  for  beta  transport  was  used  along 
with  the  ionization  and  recombination  rates  prescribed  by  them. 

A brief  review  of  ADVECT  and  the  beta  deposition  routines 
appended  to  it  is  given  in  this  report  along  with  some  typical  results 

provided  by  the  code.  This  section  concludes  with  a review  of  the  beta  i 

transport  treatment  (Knapp  and  Fischer,  1970)  and  the  results  of  a study 
of  its  relative  accuracy. 

5.1  ADVECT-Review 

A detailed  description  of  the  basic  ADVECT  computer  code  is 
found  in  Zalesak  and  Coffey  (1975).  Briefly,  the  interior  and  exterior 
surfaces  of  a nuclear  debris  cloud  after  s tabi lization  are  defined  as 
nets  of  tracer  particles  which  move  passively  in  a predetermined  velocity 
field.  At  each  time  step,  each  tracer  particle  is  moved  using  a second 
order  scheme  with  the  local  fluid  velocity  multiplied  by  the  time  step 
increment.  Thus  from  the  simplest  of  calculations,  the  motion  of  the 
entire  cloud  is  determined. 

Calculations  of  typical  stabilized  debris  clouds  using  ADVECT 
quickly  revealed  the  overwhelming  Importance  of  vertical  wind  shear, 
as  well  as  of  net  translation  of  the  cloud  by  the  winds.  While  the  models 
used  in  codes  like  WEPH  predicted  a constant  horizontal  radial  expansion 
of  the  cloud  in  time,  with  the  center  of  mass  remaining  at  the  release 
point,  ADVECT  revealed  a rapid  vertical  shearing  of  the  cloud. 
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predominantly  in  the  east-west  direction.  At  times  later  than  a few 


hours,  the  cloud  gives  the  appearance  of  a long  worm-like  structure. 

In  addition,  large  excursions  of  the  center  of  mass  of  this  structure 
from  the  initial  release  point  were  not  uncommon. 

As  an  example,  Fig.  S .1  shows  the  computed  position  and  con- 
figuration of  clouds  2k  hours  after  release  over  the  central  United 
States  (Uo°N,  265°E)  at  75  km  altitude  on  k different  days  of  the  year 
(Jan.  1,  Apr.  1,  July  1,  Oct.  1,  respectively)  at  12:00  noon.  The 
clouds  were  originally  centered  ellipsoids  20  km  thick  vertically  and 
6o  km  in  horizontal  diameter.  The  box  in  the  map  of  each  plot  shows 
the  location  of  the  longitude-latitude-altitude  "box"  enclosing  the 
cloud  given  by  ADVECT.  By  contrast,  the  circle  shows  the  cloud  location 
given  by  the  model  in  WEPH.  The  plots  below  the  map  in  each  figure 
display  views  of  the  aforementioned  longitude-latitude-altitude  box 
enclosing  the  ADVECT  cloud  from  above  and  from  the  south.  Note  that 
the  abscissa  and  ordinate  are  not  to  scale,  as  can  be  seen  from  the 
cloud-box  dimensions  given  below  the  two  plots. 

5.2  Addition  of  Beta  Deposition  Routines  to  ADVECT 

From  the  results  of  the  previous  section,  it  is  obvious  that 
one  would  expect  differences  in  communication  path  conditions  between 
ADVECT  and  the  models  used  in  WEPH.  In  order  to  quantify  these  differ- 
ences, it  was  necessary  to  add  beta  particle  deposition  and  chemistry 
routines  to  ADVECT.  We  describe  the  physics  briefly: 

For  several  days  after  the  release  a nuclear  debris  cloud  is 
an  intense  source  of  beta  particles.  These  beta  particles  are  constrained 
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Fig.  5.1  — Cloud  configurations  and  positions  24  hours  after  release  for  clouds 
released  at  75  km  altitude  over  central  United  States  on  dates  Jan.  1,  Apr.  1, 
July  1,  cuid  Oct.  1.  The  circles  on  the  map  represent  the  results  of  a horizontal 
radial  expansion  about  the  release  point,  while  the  rectangles  represent  the  results 
of  ADVECT. 


to  move  along  geomagnetic  field  lines.  Half  the  beta  particles  produced 
along  a given  field  line  are  assumed  to  escape  upward;  the  other  half 
move  downward  into  the  dense  atmosphere  below,  and  are  absorbed.  The 
absorption  process  produces  enhanced  ionization. 

At  equilibrium  the  intensity  of  the  induced  ionization  is  de- 
termined by  local  ion-electron  recombination  rates  and  the  external 
electron  deposition  rate.  The  net  result  is  a quasi-equilibrium  electron 
density  enhancement  below,  the  degree  depending  on  beta  particle  flux, 
cloud  altitude,  effective  dip  angle,  time  of  day,  latitude,  and  the 
altitude  in  question. 

The  first  step  toward  the  computation  of  electron  density 
enhancements  was  the  development  of  a subroutine  which  integrated  the 
local  beta  particle  production  rate  along  geomagnetic  field  lines 
through  the  cloud.  To  do  this  required  a more  precise  definition  of 
the  surface  of  the  cloud.  A surface  is  defined  in  our  model  as  the 
sum  of  the  subsurfaces  of  an  n x m "net"  of  tracer  points.  Thus  a 
subsurface  is  the  interior  of  a figure  formed  by  connecting  point  (i,j) 
to  i i + 1,  j)  to  ‘i  +1,  j + l1  to  (,i,  j + l)  to  vi,j),  Is  1 £ n, 

1 ^ J S m.  However,  these  four  points  do  not  define  a unique  surface. 

The  ambiguity  is  removed  by  subdividing  each  of  these  subsurfaces  with  a 
line  from  (i,j)  to  (i  + 1,  j + 1),  thus  forming  two  planar  triangular 
subdivisions  in  each  subsurface.  The  integrated  beta  particle  production 
rate,  or  beta  particle  flux,  is  found  by  computing  the  intersection,  if  any, 
of  a geomagnetic  field  line  with  each  subsurface,  multiplying  the  production 
rate  between  intersections  by  the  distance  between  intersections,  and 
summing  over  all  intersection  pairs.  The  local  production  rate  itself 
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depends  on  the  fission  yield  of  the  nuclear  device  in  question  and  on 


the  time  elapsed  after  the  burst. 

Once  the  beta  particle  flux  along  a field  line  is  known,  the 
electron  density  enhancement  at  any  point  along  that  field  line  can 
be  computed  using  the  prescriptions  given  in  Knapp  and  Fischer  (1972). 
Further  prescriptions  in  Knapp  and  Fischer  can  be  used  to  compute  the 
non-deviative  absorption  of  the  enhanced  region  due  to  electron-neutral, 
ion-neutral,  and  ion-electron  collisions  for  a given  electromagnetic 
wave  frequency. 

Figures  5-2  - show  the  results  from  ADVECT  for  a 1 megaton 

fission  yield  cloud  released  at  95  km  altitude  over  Johnston  Island 
(169.9'’  W,  16.5‘’n)  on  January  1,  6 hours  after  release.  Displayed  are 
10  MHz  wave  absorption  contours  (db/km)  as  a function  of  longitude 
(abscissa'  and  latitude  iordinate)  in  degrees  for  altitudes,  65,  75> 
and  85  km  respectively.  Figures  5*5  - 5-7  display  the  same  results  for 
a WEPH-type  code.  In  addition,  Figures  5-S  and  5*9  plot  the  latitude 
and  longitude  integrated  attenuations,  respectively,  for  both  the  ADVECT 
and  WEPH-type  code  calculations  of  the  above  problem.  These  plots  give 
the  total  signal  attenuation  in  db  for  10  MHz  signals  along  all  east- 
west  and  north-south  paths  through  the  clouds  at  each  of  the  three 
altitudes. 

It  can  be  seen  that  the  vertical  extent  of  the  significant 
signal  attenuation  is  considerable  (~30  km^  and  that  is  fact  the 
attenuation  region  presents  itself  as  a rather  thick  "curtain"  through 
which  signals  may  or  may  not  have  to  pass,  depending  on  the  exact  geometry 


of  the  signal  path.  Note  that  quite  significant  attenuation  can  take 
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Fig.  5.4  — Sii-iie  as  Fig.  5.2  but  at  altitude  85  km 
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LONGITUDE  - DEGREES 

Fig.  5.8  — Toted  10  MH^  wave  attenuation  (db)  for  north-south  directed  signal 
path  vs  longitude  of  signal  path  in  degrees.  Shown  are  plots  for  both  the  AD- 
VKCT  and  VVEPH-like  calculations  at  65,  75  and  85  km  signal  path  altitude. 
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Fig.  5.9  — Total  10  MH.^  wave  attenuation  (db)  for  east-west  directed  signal  paths 
vs.  latitude  of  signal  path  in  degrees.  Shown  are  plots  for  both  the  ADVECT  and 
WEPH-like  calculations  at  65,  75  and  85  km  signal  path  altitude. 
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place  even  at  this  late  time.  Also  note  the  dependence  of  the  total 
attenuation  on  the  path  direction  (WEPH  models  would  compute  no  path 
dependence  for  paths  through  the  cloud  center).  It  can  be  seen  that 
the  large  differences  in  total  wave  attenuation  along  the  same  signal 
path  are  due  to  two  distinct  causes:  l)  The  cloud  center  is  displaced 

in  ADVECTj  whereas  it  remains  fixed  in  the  WEPH  model;  and  2)  the  en- 
hancement regions  in  ADVECT  are  ribbon- like,  not  the  pancakes  assumed 
in  WEPH.  Thus  a signal  path  that  WEPH  determines  to  be  inoperative 
could  in  fact  be  wide  open,  and  vice-versa. 

5.3  Accuracy  of  the  Beta  Transport  Scheme 

The  prescription  used  to  determine  the  free  electron  production 
rate,  i.e.  the  analytic  expressions  for  the  beta  induced  ionization  of 
atmosphere  found  in  Knapp  and  Fischer  [19J0) , has  the  advantage  of  being 
computational ly  efficient.  We  are  interested  here  in  summarizing  the 
physical  assumptions  made  in  arriving  at  this  analytical  approximation 
and  comparing  the  results,  for  typical  cases,  produced  with  this  approxi- 
mation with  those  yielded  by  a computationally  less  tractable  method 
of  known,  high  accuracy. 

In  both  treatments  betas  are  assumed  to  be  emitted  isotropically 
and  to  follow  spiral  paths  along  geomagnetic  field  lines.  The  following 
assumptions  are  also  made  in  the  Knapp  and  Fischer  scheme: 


1) 


The  beta  energy  deposition  coefficient 
beta  energy  over  the  range  involved. 


is  independent  of 
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2)  The  beta  emission  spectrum  can  be  represented  by 

F(E)  = — ^ e - E/'  betas  / MeV  - beta 

where  E„  is  the  average  beta  energy  and  F(E)dE  represents 

D 

the  fraction  of  betas  with  the  energy  range  between  E and 
E + dE. 

3)  Energy  absorbed  from  betas  is  localized  in  the  region  of  the 
energy  absorbing  encounter. 

4)  There  is  no  change  in  the  pitch  angle  distribution  due  to 
atmospheric  collisions. 

5)  Ionization  due  to  electrons  trapped  in  the  geomagnetic  field 
is  negligible. 

Under  these  assumptions  the  ionization  rate  per  unit  volume  of 
atmosphere  can  be  expressed  as 


qg  = 


1 


F(U,  B) 
uBs  /B 


mu 


du 


(5.1) 


where  F(u,B)  du,  is  the  number  of  betas  with  pitch  angle  cosines 
between  u and  p,  + du  , n^  is  the  number  of  ion  pairs  created 
per  unit  energy  deposited,  B^  and  B are  the  geomagnetic  intensities 
at  the  source  and  deposition  points  respectively,  p is  the  atmospheric 
mass  density  at  the  deposition  point,  and  m is  the  atmospheric 
mass  penetrated  along  the  geomagnetic  field  line  between  the  source 
and  deposition  points. 
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m is  given  by: 


P - P 


g sin  cp 


where  cp  is  the  magnetic  dip  angle  and  P and  Ps  are  atmospheric 
pressures  at  the  deposition  and  source  points.  Equation  5-?  derives 
from  the  assumption  of  hydrostatic  equilibrium. 

F (a,B),  the  beta  pitch  angle  distribution  function,  can 
be  evaluated  by  noting  that  a beta  will  reflect  or  mirror  when 


s 1- 


Slnce  emission  at  the  source  is  isotropic  we,  have 


where  N is  the  number  of  betas  emitted  per  unit  time,  per  unit 
volume.  Under  the  assumption  of  negligible  pitch  angle  scattering 
from  atmospheric  collisions  the  invariance  of  beta  magnetic  moment 
leads  to 


sin  9s 
/B_ 


sin  0 


From  Eqs.  v5.5)>  (5.^''  and  (5.5) 


F(u,B) 


B 

s 

B 


/I 


u 


B 


(5.6) 


By  substituting  Expressions  (5.2)  and  (5.6'!  in  Equation  (5.1'  and  inter 
grating  over  p,  the  final  expression  for  q^  becomes 


^ Eg  /sin  cp  ^ 


(5.T) 


where  Ej^  is  the  exponential  integral  defined  by; 


E^  (x) 


e *■/!  dt 


An  Alternative  Energy  Deposition  Scheme 

The  expression  for  the  ionization  rate  due  to  the  passage  of 

betas  through  the  atmosphere  given  in  Eq.  (5.1)  was  derived  by  assuming 

a constant  energy  deposition  coefficient,  u-  , and  integrating  the 

P 

exponential  energy  spectrum  from  E = Uq  tn/u,  to  E = p-  m/u  , 

P P 

where  m is  the  total  mass  penetrated  along  a geomagnetic  field  line. 
In  general,  Uq  is  not  a constant  but  an  energy  dependent  variable 
arising  from  discrete  energy  losses  suffered  by  an  electron  in  an 
essentially  random  series  of  collisions  with  atmospheric  atoms. 

These  collisions  with  atomic  electrons  and  nuclei  result  also  in 
changes  of  direction,  a factor  which  makes  the  total  path  length 
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and  net  depth  of  penetration  unequal.  If  the  electron  energy  is 
considerably  greater  than  the  binding  energies  of  the  atomic  electrons 
encountered,  then  the  details  of  atomic  structure  can  be  ignored  and 
appropriate  averages  over  various  atomic  characteristics  may  be  used. 

The  transport  and  energy  deposition  problem  can  be  treated  within 
either  a Monte-Carlo  or  a diffusion  equation  framework,  the  latter 
being  computationally  more  efficient  if  somewhat  less  flexible. 

The  virtue  of  going  from  a constant  u-  model  to  a more 

a 

physically  realistic  one  may  be  seen  more  easily  if  we  first  consider 

the  problem  of  a planar  source  embedded  in  an  infinite  ocean  of  air, 

emitting  monoenerget ic , monodirectional  electrons.  Energy  deposition 

profiles  for  0.2,  0.4,  1.0  and  2.0  MeV  electrons  from  transport 

calculations  by  Spencer  (1959)  are  shown  in  Fig.(5.10)  along  with  constant 

2 

)j,-  profiles.  The  value  of  u,_  is  taken  as  2 MeV-cm  /gm,  the 
9 0 

optimum  value  used  in  the  Knapp-Fischer  scheme.  The  plots  show  that 
the  constant  coefficient  assumption  should  cause  an  underestimate 
in  the  deposition  from  low  energy  electrons  at  high  altitudes  and 
an  overestimate  of  the  amount  of  energy  reaching  the  lower  altitudes. 

Most  of  the  beta  energy  in  the  assumed  exponential  spectrum  falls 
within  the  0.2  - 2.0  MeV  range  of  these  plots  since  the  average 
energy,  E.  , varies  from  1.0  MeV  a few  minutes  after  detonation  to 
0.7  MeV  a few  hours  later. 

The  choice  of  a suitable  electron  transport  scheme  was  made  on 
the  basis  of  the  availability  of  an  efficient  computer  code  (Strick- 


land et . al.,  (1975)  ) for  solving  the  Fokker-Planck  equation.  The 
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form  of  the  equation  applicable  here  has  the  distribution  function 
* expressed  in  terms  of  the  variables  u,  E and  r where  a is  the 
cosine  of  the  orbit  pitch  angle,  E is  the  electron  kinetic  energy 
and  is  an  effective  column  density  given  by 


dT  = n(z)  CJ^  (E)  dz 


where  n (z)  is  the  mass  density  of  the  atmosphere  and  is  the 

sum  of  the  ionization,  excitation  and  elastic  collision  cross  sections. 
The  equation  is  expressed  in  the  form: 


(t  , E , jj. ) 


q(E). 

2o  (E) 


i c 


n 2,  M 1 


1 ^ 
(T  (E)  dE 


( L (E)  $ ) 


where  Z (E)  and  L (E)  are  the  momentum  transfer  cross  section  and 
energy  loss  functions  respectively. 

The  original  low  energy  version  of  the  Fokker-Planck  solver  was 
modified  by  the  substitution  of  relativistic  cross  sections  for  non- 
relativistic  ones.  The  loss  of  energy  suffered  by  an  electron  per 
unit  path  length  as  a result  of  collisions  with  orbital  electrons 
of  the  medium  is  given  by  the  formulation  of  Rohrlich  and  Carlson 
(1954)  of  the  Bethe-Block  theory: 
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where  L (T)  is  the  average  energy  loss  per  unit  path  length 

2 2 
(MeV-cm  /gm)  ; E is  the  electron  kinetic  energy  (MeV)  , me  is  the 

2 

electron  rest  mass  energy,  T = E/mc  ; B = v/c  and  , and 

are  the  atomic  weight,  atomic  number,  and  the  fraction  by  mass  of 
the  i'th  constituent  of  the  medium.  I is  the  mean  excitation  energy 
of  the  medium,  taken  as  96.8  eV  for  air. 

The  momeiitum  transfer  cross  section  was  calculated  using  the 
screened  Rutherford  cross  section  of  Goudsmit  and  Saunderson  (1940) 
multiplied  by  a spin  factor  evaluated  by  McKinley  and  Feschback  (1948). 
This  differential  cross  section  for  scattering  into  solid  angle 
2tt  sin  0 d0  is  given  by: 


dcT  _ 

z e 

1 

i. 

dn 

2 4 

m V 

(1-cos  0 + 

2T])^ 

TTQfB 

72 

(1-cos  0 )^ 

(1  - 

(1-cos  9 ) 


(1-cos  0)^  , 
/2  ^ 


I 
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The  code  incorporating  these  modifications  was  checked  by 
duplicating  a case  studied  by  Berger,  Seltzer,  and  Maeda  (1970) 
involving  energy  deposition  in  the  atmosphere  by  electrons  with  an 
exponential  energy  spectrum  with  average  energy  200  KeV.  The  Berger 
et . al.  study  utilized  a Monte-Carlo  code.  The  deposition  profile 
provided  by  the  Fokker-Planck  code  duplicated  the  Monte-Carlo  results 
within  the  bounds  of  the  accuracy  of  the  graphical  display  provided 
by  Berger  et,  al. 

Conclusions 

In  order  to  access  the  limitations  of  the  simple  constant 
energy  deposition  scheme  it  was  compared  to  Fokker-Planck  results 
for  two  cases:  Eg  = 0.8  MeV  and  Eg  = 1.0  MeV.  The  first  is 

characteristic  of  a beta  spectrum  about  an  hour  after  detonation, 
the  second  about  a few  seconds  after  detonation.  The  energy  deposition 
results,  plotted  in  Figs. (5.11/  and  (5.12'  show  that  the  simple  model  agrees 
quite  well  with  the  Fokker-Planck  curves  between  60  km  and  90  km 
altitude,  with  the  disparity  typically  30%  or  smaller.  At  50  km 
and  below,  where  the  deposition  falls  off,  this  disparity  between 
the  treatments  increases  rapidly  to  an  order  of  magnitude  and 
greater,  reflecting  differences  seen  in  the  deposition  curves  for 
monochromatic  betas.  These  results  are  for  an  exponential  spectrum, 
but  should  hold  equally  for  an  empirical  beta  spectrum. 
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Section  6 


SUMMARY 

Radioactive  debris  clouds  from  high  altitude  nuclear  explosions 
ionize  air  in  the  magnetically  conjugate  regions  at  lower  altitudes  pri- 
marily through  the  emission  of  energetic  betas  which  undergo  ionizing 
collisions  with  atmospheric  constituents.  The  ionized  regions  thus 
formed  constitute  a source  of  interference  for  satellite  communications. 
Since  the  betas  are  constrained  to  spiral  along  geomagnetic  field  lines^ 
the  beta  flux  pattern,  to  a first  approximation,  has  the  same  shape  as 
the  debris  cloud.  The  debris  cloud,  on  a time  scale  of  hours,  changes 
shape  and  location  through  advection  by  mesospheric  winds.  Past  NRL 
results  indicate  that  short  scale  length  shears  observed  in  the  meso- 
spheric circulation  typically  give  rise  to  highly  elongated  cloud  con- 
figurations within  hours  after  burst.  The  understanding  of  such  late 
time  HANE  effects  is,  therefore,  critically  dependent  on  knowledge  of 
upper  atmosphere  wind  patterns. 

This  past  year,  NRL  has  completed  an  analysis  of  observational 
data  on  mesospheric  winds,  in  particular,  the  climatological  model  given 
by  Groves,  CIRA  il972).  The  modelvas  found  to  have  large  systematic 
errors  attributable,  primarily,  to  tidal  components  and  planetary  waves 
which  remain  unresolved  in  the  sparse  data  field.  NRL  has  developed 
theoretical  dynamical  models  of  the  upper  atmosphere  which,  in  conjunction 
with  observational  data,  will  provide  a reliable  predictive  capability 
for  the  mesospheric  circulation. 

The  circulation  of  the  upper  atmosphere  is  thermally  driven 
by  Insolation  which  varies  in  time  and  space.  Radiative  heating  functions 
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have  been  developed  to  replace  earlier,  less  accurate  ones  by  Leovy 
used  in  the  NRL  linear  model  of  the  mean  mesospheric  circulation.  The 
heating  functions  are  produced  by  a computer  program  which  can  be  used 
to  provide  either  mean  seasonal  heating  rates  for  use  in  climatological 
forecasts  or  point-by-point  local  heating  rates  suitable  for  solar  tidal 
calculations  (diurnal  and  semi-diurnal  time  scales).  The  new  heating 
function  has  been  incorporated  into  the  NRL  linear  mesospheric  model 
and  test  runs  have  indicated  that  the  model  accurately  simulates  the 
observed  mean  zonal  flow  for  both  solstitial  and  equinoctial  conditions. 

Superimposed  on  the  mean  mesospheric  circulation  are  planetary 
waves  on  a diurnal  time  scale  which  react  non-linearly  with  each  other 
and  with  the  mean  wind  current.  This  solar-driven  tidal  fluctuation  has 
been  modeled  at  NRL  with  a non-linear  semi-spectral  numerical  scheme 
which  has  successfully  simulated  both  the  diurnal  and  semi-diurnal 
modes  of  the  tide.  This  program  has  been  improved  in  the  past  year  with 
major  changes  in  the  lower  boundary  condition  and  increased  resolution 
in  the  computational  grid.  Comparison  of  results  from  this  model  with 
analytic  results  from  Lindzen  (1971)  reveal  excellent  agreement. 

The  overall  goal  of  the  NRL  effort  is  the  development  of  a 
numerical  model  which,  given  wind  and  bomb  data,  tracks  the  advection 
of  debris,  computes  the  equilibrium  distribution  of  beta  induced  ioni- 
zation at  specified  time,  and  then  determines  coimtunicat ion  signal 
attenuation  by  ray  tracing. This  program  has  been  completed  with  the 
incorporation  of  beta  deposition  routines  based  on  analytic  expressions 
given  by  Knapp  and  Fischer  (1970).  iTiese  approximate  expressions  have 
been  found  to  agree  within  + with  a more  rigorous  NRL  electron 
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transport  code  based  on  the  Fokker-Planck  equation,  in  the  regions 
of  peak  ionization. 
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OICY  ATTN  W14  PAT  CLARK 
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OJCS/ J-5 
THE  PENTAGON 
WASHINGTON,  O.C.  20301 
(OPERATIONS) 

OICY  ATTN  WWMCCS  EVAL  OFC  MR.  TOMA 


ASD  CC51)  (SYSTEMS) 

3D224,  THE  PENTAGON 
WASHINGTON,  D.C.  20301 

OICY  ATTN  OR.  M.  EPSTEIN 
OICY  DR.  J.  BABCOC< 


WWMCCS  SYSTEM  ENGINEERING  ORG 
WASHINGTON,  O.C.  20305 

OICY  ATTN  R.  L.  CRAWFORD 


commander/director 

ATMOSPHERIC  SCIENCES  LABORATORY 
U.S.  ARMY  electronics  COMMAND 
WHITE  SANDS  MISSILE  RANGE,  NM  88002 

OICY  ATTN  DRSEL-BL-SY-S  F.  E.  NILES 


DIRECTOR 

BMD  ADVANCED  TECH  CTR 
HUNTSVILLE  OFFICE 
P.  0.  BOX  1500 
HUNTSVILLE,  AL  35807 
OICY  ATTN  ATC-T 
OICY  ATTN  ATC-0 
OICY  ATTN  ATC-R 


MELVIN  T.  CAPPS 
W.  DAVIES 
DON  RUSS 


PROGRAM  MANAGER 
BMD  PROGRAM  OFFICE 
5001  EISENHOWER  AVENUE 
ALEXANDRIA,  VA  22333 

OICY  ATTN  DACS-BMT  JOHN  SHEA 


CHIEF  C-E  SERVICES  DIVISION 
U.S.  ARMY  COMMUNICATIONS  CMO 
PENTAGON  RM  1B259 
WASHINGTON,  D.C.  20310 

OICY  ATTN  CC-OPS-CE 


-3- 


COMMANDER 

HARRY  DIAMOND  LABORATORIES 
2800  POWDER  MILL  ROAD 
ADELPHI,  MD  20783 

(CNWDI-INNER  ENVELOPE:  ATTN:  DRXDO-RBH) 

OICY  ATTN  MILDRED  H.  WEINER  DRXDO-TI 
OICY  ATTN  DRXDO-RB  ROBERT  WIlLIAMS 
OICY  ATTN  DRXDO-NP  FRANCIS  N.  WIMENITZ 
OICY  ATTN  DRXDO-NP  CYRUS  MOAZEO 


DIRECTOR 

TRASANA 

WHITE  SANDS  MISSILE  RANGE,  NM  88002 
OICY  ATTN  ATAA-SA 
OICY  ATTN  TCC/F.  PAYAN  JR 
OICY  ATTN  ATAA-TAC  LTC  JOHN  HESSE 


COMMANDER 

U.S.  ARMY  COMM-ELEC  ENGRG  INSTAL  AGY 
FT.  HUACHUCA,  AZ  85615 

OICY  A-"N  EED-PED  GEORGE  LANE 


COMMANDER 

U.S.  ARMY  ELECTRONICS  COMMAND 
FORT  MONMOUTH,  NJ  07703 

OICY  ATTN  DRSEL-NL-RD  H.  S.  BENNET 
OICY  ATTN  DRSEL-PL-ENV  HANS  A.  BOMKE 


COMMANDER 

U.S.  ARMY  FOREIGN  SCIENCE  S TECH  CTR 
220  7TH  STREET,  NE 
CHARLOTTESVILLE,  VA  22901 

OICY  ATTN  P.  A.  CROWLEY 
OICY  ATTN  R.  JONES 


COMMANDER 

U.S.  ARMY  MATERIEL  DEV  6 READINESS  CMD 
5001  EISENHOWER  AVENUE 
ALEXANDRIA,  VA  22333 

OICY  ATTN  DRCLOC  J.  A.  BENDER 
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COMMANDER 

U.S.  ARMY  NUCLEAR  AGENCY 
FORT  3LIS5,  TX  79916 

OICY  ATTN  MONA-WE  J.  3ER8ERET 


DIRECTOR 

U.S.  ARMY  ballistic  RESEARCH  LABS 
ABERDEEN  PROVING  GROUND,  MD  21005 
OICY  ATTN  LAWRENCE  J.  PUCKETT 


commander 

U.S.  ARMY  SATCOM  AGENCY 
FT.  MONMOUTH,  NJ  07703 

OICY  ATTN  DOCUMENT  CONTROL 


commander 

U.S.  ARMY  MISSILE  INTELLIGENCE  AGENCY 
REDSTONE  ARSENAL,  AL  35809 
OICY  ATTN  JIM  GAMBLE 


CHIEF  OF  NAVAL  OPERATIONS 
NAVY  DEPARTMENT 
WASHINGTON,  D.C.  20350 

OICY  ATTN  OP  943  LCDR  HUFF 
OICY  ATTN  ALEXANDER  BRANDT 
OICY  ATTN  RONALD  E.  PITKIN 


CHIEF  OF  NAVAL  RESEARCH 
NAVY  DEPARTMENT 
ARLINGTON,  VA  22217 

OICY  ATTN  CODE  418 
OICY  ATTN  CODE  461 


COMMANDER 

naval  ELECTRONIC  SYSTEMS  COMMAND 
NAVAL  ELECTRONIC  SYSTEMS  CMD  HQS 
WASHINGTON,  D.C.  20360 

OICY  ATTN  NAVALEX  034  T BARRY  HUGHES 

OICY  ATTN  PME  117-T  SATELLITE  COMM  PROJECT  OFF 

OICY  ATTN  PME  117 


COMMANDING  OFFICER 
NAVAL  intelligence  SUPPORT  CTR 
4301  SUITLAND  road,  bldg.  5 
WASHINGTON,  D.C.  20390 

OICY  ATTN  MR.  DUBBIN  STIC  12 
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COMMANDER 


naval  ocean 

SYSTEMS  CENTER 

SAN  DIEGO, 

CA  92152 

0 ICY 

ATTN  W I LL 1 AM  F . 

MOLER 

0 ICY 

ATTN  CODE  0230 

C.  BAGGETT 

03CY 

ATTN  CODE  2200 

0 ICY 

ATTN  R.  EASTMAN 

D [ RECTOR 

NAVAL  RESEARCH  LABORATORY 


WASHINGTON 

, D.C  , 

. 2J375 

0 ICY 

ATTN 

CODE 

7 7 0 0 

TIMOTHY  P.  COFFEY 

0 ICY 

ATTN 

CODE 

5460 

ELECTRCMAG  PROP  BR 

OICY 

ATTN 

CODE 

5430 

03CY 

ATTN 

CODE 

770  1 

JACK  D.  BROWN 

0 ICY 

ATTN 

HDQ 

COMM 

DIR  BRUCE  WALD  (CODE 

OICY 

ATTN 

CODE 

5465 

PROP  applications 

0 ICY 

ATTN 

CODE 

546  1 

TRANS  lONO  PROP 

OICY 

ATTN 

CODE 

7 7 50 

COMMANDER 

NAVAL  SPACE  SURVEILLANCE 
OAHLGREN,  VA  22448 

SYSTEM 

OICY 

ATTN 

CAPT 

J.  H 

. BURTON 

5400) 


COMMANDER 

NAVAL  SURFACE  WEAPONS  CENTER 
WHITE  OAK,  SILVER  SPRING,  MD  20910 

OICY  ATTN  CODE  WA501  NAVY  NUC  PRGMS  OFF 


DIRECTOR 

STRATEGIC  SYSTEMS  PROJECT  OFFICE 
NAVY  DEPARTMENT 
WASHINGTON,  O.C.  20376 
OICY  ATTN  NSP-2141 

OICY  ATTN  NSSP-2722  FRED  WIMBERLY 


NAVAL  SPACE  SYSTEM  ACTIVITY 
P.  0.  BOX  92960 
WORLOWAY  POSTAL  CENTER 
LOS  ANGELES,  CALIF.  90009 

OICY  ATTN  A.  B.  HAZZARD 


COMMANDER 

ADC/DC 

ENT  AFB,  CO  8 '912 

OICY  ATTN  DC  MR.  LONG 
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COMMANDER 
AOCOM/XPD 
ENT  AFB,  CO  8C912 

OICY  ATTN  XPQDQ 


AF  GEOPHYSICS  LABORATORY.  AFSC 


MANSCOM  AFB,  MA  017 
OICY  ATTN  OPR 
OICY  ATTN  OPR 
OICY  ATTN  LKB 
OICY  ATTN  OPR 
OICY  ATTN  SUOL 
OICY  ATTN  PHP 
OICY  ATTN  PHD 
OICY  ATTN  PHD 


1 

HAROLD  GARDNER 
JAMES  C.  ULWICK 
KENNETH  S.  W.  CHAMPION 
ALVA  T.  STAIR 
RSCH  LIB 
JULES  AARONS 
JURGEN  BUCHAU 
JOHN  P.  MULLEN 


AF  WEAPONS  LABORATORY,  AFSC 
KIRTLAND  AFB,  NM  87117 


AFTAC 


0 

ICY 

ATTN 

SUL 

0 

ICY 

ATTN 

CA 

ARTHUR 

H 

. GUENTHEl 

0 

ICY 

ATTN 

DYC 

CAPT  L 

, 

WITTWER 

0 

ICY 

ATTN 

SAS 

JOHN 

M . 

KAMM 

0 

ICY 

ATTN 

DYC 

CAPT 

MARK  A.  FRY 

CK  AFB 

, PL 

529  25 

0 

ICY 

ATTN 

TF/MAJ  WIL 

EY 

0 

ICY 

ATTN 

TN 

AIR  FORCE  AVIONICS  LABORATORY,  AFSC 
WRIGHT-PATTERSON  AFB,  OH  95433 
OICY  ATTN  AAD  WADE  HUNT 
OICY  AiTN  AFAL  AA3  H.  M.  HARTMAN 
OICY  ATTN  AAD  ALLEN  JOHNSON 


HEADQUARTERS 

ELECTRONIC  SYSTEMS  DIVI5ION/XR 
HANSCOM  AFB,  MA  01731 

OICY  ATTN  XRC  LTC  J.  MORIN 
OICY  ATTN  XRE  LT  MICHAELS 


HEADQUARTERS 

ELECTRONIC  SYSTEMS  OlVlSiON/YS 
HANSCOM  AFB,  MA  01731 
OICY  ATTN  'YSEV 


r 


COMMANDER 
FOREIGN  TE 
WRIGHT-PAT 
0 ICY 
0 ICY 


CHNOLOGY  DIVISION,  AFSC 
TER50N  AFB,  OH  45433 
ATTN  NICD  LIBRARY 
ATTN  ETD  B.  L.  BALLARD 


HQ  USAF/RD 

WASHINGTON,  D.C.  20330 
OICY  ATTN  RDQ 


COMMANDER 

ROME  AIR  DEVELOPMENT  CENTER,  AFSC 
GRIFFISS  AFB,  NY  13440 

OICY  ATTN  EMTLD  DOC  LIBRARY 
OICY  ATTN  OCSE  V.  COYNE 


SAMSO/SZ 

POST  OFFICE  BOX  92960 
WORlDWAY  POSTAL  CENTER 
LOS  ANGELES,  CA  90009 
(SPACE  DEFENSE  SYSTEMS) 

OICY  ATTN  SZJ  MAJOR  LAWRENCE  DOAN 


COMMANDER  IN  CHIEF 
STRATEGIC  AIR  COMMAND 
OFFUTT  AFB,  NB  68113 

OICY  ATTN  XPFS  MAJ  BRIAN  G.  STEPHAN 
OICY  ATTN  ADWATE  CAPT  BRUCE  BAUER 
OICY  ATTN  NRT 


HEADQUARTERS 

ELECTRONIC  SYSTEMS  DIVISION  (AFSC) 
HANSCOM  AFB,  MA  01731 

OICY  ATTN  JIM  DEAS 


SAMSO/ YA 
P.  0.  BOX  92960 
worldway  postal  CENTER 
LOS  ANGELES,  CA  90009 

OICY  ATTN  YAT  CAPT  L.  8LACKWELDER 


SAMSO/SK 
P.  0.  BOX  92960 
WORLDWAY  postal  CENTER 
lOS  ANGELES,  CA  90009 

OICY  ATTN  SKA  LT  MARIA  A.  CLAVlN 


-8- 


SAMSO/MN 

NORTON  AFB,  CA  92419 
(MINUTEMAN) 

OICY  ATTN  MNNL  LTC  KENNEDY 


COMMANDER 

ROME  AIR  DEVELOPMENT  CENTER,  AFSC 
HANSCOM  AFB,  MA  01753 

OICY  ATTN  ETEI  A.  LORENTZEN 
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U.S.  ENERGY  RSCH  AND  DEV  ADMIN 


EGSG,  INC. 

LOS  alamos  DIVISION 

P.  0.  BOX  809 

LOS  alamos,  NM  85544 

OICY  ATTN  JAMES  R.  BREEDLOVE 


university  of  CALIFORNIA 
LAWRENCE  LIVERMORE  LABORATORY 
P.  0.  BOX  808 
LIVERMORE,  CA  94550 

OICY  ATTN  TECH  INFO  DEPT  L-3 
OICY  ATTN  RONALD  L.  OTT  L-531 
OICY  ATTN  DONALD  R.  DUNN  L-156 
OICY  ATTN  RASPH  S.  HAGER  L-31 


LOS  ALAMOS  SCIENTIFIC  LABORATORY 
P.  0,  BOX  1663 


ALAMOS 

, NM 

8 7 545 

OICY 

ATTN 

DOC 

CON 

FOR 

ERIC 

LINDMAN 

0 ICY 

ATTN 

DOC 

CON 

FOR 

R . F , 

. TASCHEK 

0 ICY 

ATTN 

DOC 

CON 

FOR 

ERIC 

JONES 

0 ICY 

ATTN 

DOC 

CON 

FOR 

JOHN 

S.  MALIK 

0 ICY 

ATTN 

DOC 

CON 

FOR 

MARTIN  TIERNEY  J-IO 

OICY 

ATTN 

DOC 

CON 

FOR 

JOHN 

Z INN 

SANDIA  LABORATORIES 
P.  0.  BOX  5800 
ALBUQUERQUE,  NM  87115 

OICY  ATTN  DOC  CON  FOR  J.  P.  MARTIN  ORG  1732 

OICY  ATTN  DOC  CON  FOR  W.  D.  BROWN  ORG  1353 

OICY  ATTN  DOC  CON  FOR  A.  DEAN  THORNBROUGH  ORG  1245 

OICY  ATTN  DOC  CON  FOR  T.  WRIGHT 

OICY  ATTN  DOC  CON  FOR  D.  A.  DAHLGREN  ORG  1722 
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other  government 


DEPARTMENT  OF  COMMERCE 
OFFICE  OF  TELECOMMUNICATIONS 
INSTITUTE  FOR  TELECOM  SCIENCE 
BOULDER,  CO  80302 

OICY  ATTN  WILLIAM  F.  UTLAUT 
OICY  ATTN  G.  REED 


N :iONAL  OCEANIC  & ATMOSPHERIC  AOMIN 
ENVIRONMENTAL  RESEARCH  LABORATORIES 
DEPARTMENT  OF  COMMERCE 
BOULDER,  CO  80302 

OICY  ATTN  JOSEPH  H.  POPE 
OICY  ATTN  RICHARD  GRUBB 
0 ICY  ATTN  C . L . RUFFNACH 


NASA 

GODDARD  SPACE  FLIGHT  CENTER 
GREENBELT,  MD  20771 

OICY  ATTN  ATS-6  OFC  P.  CORRIGAN 


; 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS 


AEROSPACE  CORPORATION 

P.  0.  BOX  92957 

LOS  ANGELES,  CA  90009 


OICY 

ATTN 

IRVING 

M.  GARFUNKEL 

0 ICY 

ATTN 

T . 

M . 

SALMI 

0 ICY 

ATTN 

V . 

JOSEPHSON 

OICY 

ATTN 

s . 

P . 

BOWER 

OICY 

ATTN 

N . 

D. 

STOCKWELL 

0 ICY 

ATTN 

P . 

P . 

OLSEN  120  RM  2224E 

0 ICY 

ATTN 

F . 

E . 

BOND  RM  5003 

OICY 

ATTN 

J . 

E . 

CARTER  120  RM  2209 

OICY 

ATTN 

F . 

A . 

MORSE  A6  RM  2407 

analytical 

SYSTEMS 

ENGINEERING  CORP 

5 OLD  CONCORD  ROAD 

BURLINGTON 

, MA 

01803 

0 ICY 

ATTN 

RADIO 

SCIENCES 

BOEING  COMPANY,  THE 
P.  0.  BOX  3707 
SEATTLE,  WA  98124 

OICY  ATTN  GLEN  KEISTER 
OICY  ATTN  D.  MURRAY 


CALIFORNIA  AT  SAN  DIEGO,  UNIV  OF 
3175  MIRAMAR  ROAD 
LA  JOLLA,  CA  92037 

OICY  ATTN  HENRY  G.  BOOKER 


BROWN  ENGINEERING  COMPANY,  INC. 
CUMMINGS  RESEARCH  PARK 
HUNTSVILLE,  AL  35807 

OICY  ATTN  ROMEO  A.  OELIBERIS 


CHARLES  STARK  DRAPER  LABORATORY,  INC. 
555  TECHNOLOGY  SQUARE 
CAMBRIDGE,  MA  02139 

OICY  ATTN  D.  B.  COX 

OICY  ATTN  J.  P.  GILMORE  MS  63 
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1 


COMPUTER  SCIENCES  CORPORATION 
P.  0.  BOX  530 
6565  ARLINGTON  8LVD 
falls  church,  VA  22046 
OICY  ATTN  H.  BLANK 
OICY  ATTN  JOHN  SPOOR 


COMSAT  LABORATORIES 
LINTHICUM  ROAD 
CLARKSBURG,  MD  20734 

OICY  ATTN  R.  R.  TAUR 


CORNELL  UNIVERSITY 

DEPARTMENT  OF  ELECTRICAL  ENGINEERING 
ITHACA,  NY  14850 

OICY  ATTN  0.  T.  FARLEY  JR 


ESL  INC. 

495  JAVA  lRIVE 
SUNNYVALE,  CA  94086 

OICY  ATTN  R.  K.  STEVENS 

OICY  ATTN  J.  ROBERTS 

OICY  ATTN  JAMES  MARSHALL 

OICY  ATTN  V.  L.  MOWER 

OICY  ATTN  C.  W.  PRETTIE 


FORD  AEROSPACE  S COMMUNICATIONS  CORP 
3939  FABIAN  WAY 
PALO  ALTO,  CA  94303 

OICY  ATTN  J.  T.  MATTINGLEY  MS  X22 


general  electric  COMPANY 
SPACE  DIVISION 
VALLEY  FORGE  SPACE  CENTER 
GODDARD  BLVD  KING  OF  PRUSSIA 
P.  0.  BOX  8555 
PHILADELPHIA,  PA  19101 

OICY  ATTN  M.  H.  BORTNER  SPACE  SCI  LAB 
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general  electric  company 

TEMPO-CENTER  FOR  ADVANCED  STUDIES 
816  STATE  STREET  (P.O.  DRAWER  QQ) 
SANTA  BARBARA,  CA  93102 
OICY  ATTN  DASIAC 
OICY  ATTN  DON  CHANDLER 
OICY  ATTN  TOM  BARRETT 
OICY  ATTN  TIM  STEPHENS 
OICY  ATTN  WARREN  S.  KNAPP 
OICY  ATTN  WILLIAM  MCNAMERA 
OICY  ATTN  B.  GAMBILL 
OICY  ATTN  MACK  STANTON 


general  research  corporation 
P.  0.  BOX  3587 
SANTA  BARBARA,  CA  93105 
i OICY  ATTN  JOHN  I SE  JR 

OICY  ATTN  JOEL  GARBARINO 


GEOPHYSICAL  INSTITUTE 
university  of  ALASKA 

FAIRBANKS,  AK  99701 

call  class  ATTN:  security  OFFICER) 

OICY  ATTN  T.  N.  DAVIS  CUNCL  ONLY) 

OICY  ATTN  NEAL  BROWN  CUNCL  ONLY) 

OICY  ATTN  TECHNICAL  LIBRARY 


GTE  SYLVANIA,  INC. 

ELECTRONICS  SYSTEMS  GRP-EASTERN  DIV 
77  A STREET 
NEEDHAM,  MA  02199 

OICY  ATTN  MARSHAL  CROSS 


INSTITUTE  FOR  DEFENSE  ANALYSES 
400  ARMY-NAVY  DRIVE 
ARLINGTON,  VA  22202 

OICY  ATTN  d.  M.  AEIN 
OICY  ATTN  ERNEST  SAUER 
OICY  ATTN  HANS  WOLFHARD 
OICY  ATTN  JOEL  8ENGSTEN 


HARRIS  CORP,  E.S.D. 

P.  0.  BOX  37 
MELBOURNE,  FL  32901 

OICY  ATTN  ADV  PROG  DEPT  DR.  CARL  DAVIS 
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A 


HSS,  INC. 

2 ALFRED  CIRCLE 
BEDFORD,  MA  01730 

OICY  ATTN  DONALD  HANSEN 


INTL  tel  S TELEGRAPH  CORPORATION 
500  WASHINGTON  AVENUE 
NUTLEY,  NU  07110 

OICY  ATTN  technical  LIBRARY 


JAYCOR 

1401  CAMINO  del  MAR 
DEL  MAR,  CA  92014 

OICY  ATTN  S.  R.  GOLDMAN 


I JOHNS  HOPKINS  UNIVERSITY 

applied  PHYSICS  LABORATORY 
• JOHNS  HOPKINS  ROAD 

I LAUREL,  MO  20810 

I OICY  ATTN  DOCUMENT  LIBRARIAN 

OICY  ATTN  THOMAS  POTEMRA 
OICY  ATTN  JOHN  DASSOULAS 


LOCKHEED  MISSILES  S SPACE  CO  INC 
P.  0.  BOX  504 
SUNNYVALE,  CA  94088 

OICY  ATTN  DEPT  60-12 
OICY  ATTN  D.  R.  CHURCHILL 


LOCKHEED  MISSILES  AND  SPACE  CO  INC 
3251  HANOVER  STREET 
PALO  alto,  CA  94304 

OICY  ATTN  MARTIN  WALT  DEPT  52-10 
OICY  ATTN  RICHARD  G.  JOHNSON  DEPT 
OICY  ATTN  BILLY  M.  MCCORMAC  DEPT  5 


KAMAN  SCIENCES  CORP 

P.  0.  BOX  7463 

COLORADO  SPRINGS,  CO  80953 

OICY  ATTN  8.  J.  BITTNER 


LINKABIT  CORP 
I*  1 045  3 ROSELLE 

SAN  DIEGO,  CA  92121 

OICY  ATTN  IRWIN  JACOBS 


2-12 

-54 
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V.l.T.  LINCOLN  LA80RA70RY 
P.  0.  BOX  73 
LEXINGTON,  MA  02173 

OICY  ATTN  LIB  A-082  FOR  DAVID  M. 
OICY  ATTN  MR.  WALDEN  X113 
OICY  ATTN  JAMES  H.  PANNELL  L-246 
OICY  ATTN  D.  CLARK 


MCDONNELL  DOUGLAS  CORPORATION 
5301  30LSA  AVENUE 
HUNTINGTON  BEACH,  CA  92647 
OICY  ATTN  N.  HARRIS 
OICY  ATTN  J.  MOULE 
OICY  ATTN  GEORGE  MROZ 
OICY  ATTN  BILL  OLSON 


MISSION  RESEARCH  CORPORATION 

735  STATE  STREET 

SANTA  BARBARA,  CA  93101 


0 ICY 

ATTN 

P.  FISCHER 

0 ICY 

ATTN 

W.  F.  CREVIER 

OICY 

ATTN 

STEVEN  L.  GUTSCHE 

0 ICY 

ATTN 

D.  SAPPENFIELD 

OICY 

ATTN 

R.  BOGUSCH 

0 ICY 

ATTN 

R.  HENDRICK 

0 ICY 

ATTN 

RALPH  KILB 

OICY 

ATTN 

DAVE  SOWLE 

0 ICY 

ATTN 

F.  FAJEN 

0 ICY 

ATTN 

M.  SCHEIBE 

0 ICY 

ATTN 

CONRAD  L.  LONGMIRE 

OICY 

ATTN 

WARREN  A.  SCHLUETE 

ITRE  CORPORATION, 

THE 

. 0.  BOX 

208 

EDFORD, 

MA  01730 

0 ICY 

ATTN  J. 

C.  KEENAN 

0 ICY 

ATTN  G. 

HARDING 

0 ICY 

ATTN  CHIEF  SCIENTIST  W 

0 ICY 

ATTN  S. 

A MORIN 

OICY 

ATTN  C . 

E.  CALLAHAN 

SEN 


PACIFIC-SIERRA  RESEARCH  CORP 
1456  CLOVERFIELD  BLVD . 

SANTA  MONICA,  CA  90404 

OICY  ATTN  F.  C.  FIELD  JR 


TOWLE 


-16- 


PHOTOMETRICS,  INC. 

442  MARRETT  ROAD 
LEXINGTON,  MA  02173 

OICY  ATTN  IRVING  L.  KOFSKY 


physical  DYNAMICS  INC. 

P.  0.  BOX  21589 
SEATTLE,  WA  98111 

OICY  ATTN  E.  J.  FREMOUW 


PHYSICAL  DYNAMICS  INC. 

P.  0.  BOX  1069 
BERKELEY,  CA  94701 

OICY  ATTN  A.  THOMPSON 
OICY  ATTN  JOSEPH  B.  WORKMAN 


R & D ASSOCIATES 

P.  0.  BOX  9695 

MARINA  DEL  REY,  CA  90291 

OICY  ATTN  RICHARD  LATTER 
OICY  ATTN  FORREST  GILMORE 
OICY  ATTN  BRYAN  GABBARD 
OICY  ATTN  william  B.  WRIGHT  JR 
OICY  ATTN  ROBERT  F.  LELEVIER 
OICY  ATTN  william  J.  KAR2AS 


RAND  CORPORATION,  THE 
1700  MAIN  STREET 
SANTA  MONICA,  CA  90406 

OICY  ATTN  CULLEN  CRAIN 
OICY  ATTN  ED  BEDROZIAN 


SCIENCE  applications,  INC. 
P.  0.  BOX  2351 
LA  JOLLA,  CA  92038 


OICY 

ATTN 

LEWIS  M.  LINSON 

0 ICY 

ATTN 

DANIEL  A.  HAMLIN 

OICY 

ATTN 

D.  SACHS 

0 ICY 

ATTN 

E.  A.  STRAKER 

0 ICY 

ATTN 

CURTIS  A.  SMITH 

0 ICY 

ATTN 

JACK  MCDOUGAL 

RAYTHEON  CO. 

528  BOSTON  POST  ROAD 
SUDBURY,  MA  01776 

OICY  ATTN  BARBARA  ADAMS 


-17- 


r 


SCIENCE  applications,  INC. 
HUNTSVILLE  DIVISION 
2109  W.  CLINTON  AVENUE 
SUITE  700 

HUNTSVILLE,  AL  35805 

OICY  ATTN  DALE  H.  DIVIS 


SCIENCE  APPLICATIONS,  INCORPORATED 
8400  WESTPARK  DRIVE 
MCLEAN,  VA  22101 

OICY  ATTN  B.  ADAMS 


STANFORD  RESEARCH  INSTITUTE 
333  RAVENSWOOD  AVENUE 
MENLO 


PARK 

, CA 

94025 

OICY 

ATTN 

DONALD  NEILSON 

OICY 

ATTN 

ALAN  BURNS 

0 ICY 

ATTN 

G.  SMITH 

0 ICY 

ATTN 

L.  L.  COBB 

OICY 

ATTN 

DAVID  A.  JOHNSON 

0 ICY 

ATTN 

WALTER  G.  CHESNUT 

0 ICY 

ATTN 

CHARLES  L.  RING 

OICY 

ATTN 

WALTER  JAYE 

OICY 

ATTN 

M.  BARON 

0 ICY 

ATTN 

RAY  L.  LEADABRAND 

SYSTEM  DEVELOPMENT  CORPORATION 
4130  LINDEN  AVENUE 
DAYTON,  OH  45432 

OICY  ATTN  F.  G.  MEYER 


technology  INTERNATIONAL  CORP 
75  WIGGINS  AVENUE 
BEDFORD,  MA  01730 

OICY  ATTN  W.  P.  BOQUIST 


TRW  DEFENSE  6 SPACE  SYS  GROUP 

ONE  SPACE  PARK 

REDONDO  BEACH,  CA  90278 

OICY  ATTN  R.  K.  PLEBUCH  Rl-2078 
OICY  ATTN  ROBERT  M.  WEBB  Rl-1150 


VISIOYNE,  INC. 

19  third  avenue 

NORTH  WEST  INDUSTRIAL  PARK 
BURLINGTON,  MA  01803 

OICY  ATTN  CHARLES  HUMPHREY 
OICY  ATTN  J.  W.  CARPENTER 
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